NAVAL POSTGRADUATE SCHOOL 
Monterey , California 




THESIS 



NUMERICAL STUDIES OF LOCALIZED VIBRATING STRUCTURES 
IN NONLINEAR LATTICES 

by 

Brian Russell Galvin 
March 1991 

Thesis Advisor: Bruce Denardo 

Co-Advisor Andres Larraza 



Approved for public release; distribution unlimited 




Unclassified 

SECURITY CLASSIFICATION OF THIS PAGE 



REPORT DOCUMENTATION PAGE 



Form Approved 
OMB No 0704-0188 



la REPORT SECURITY CLASSIFICATION 

UNCLASSIFIED 



lb RESTRICTIVE MARKINGS 



2a SECURITY CLASSIFICATION AUTHORITY 



2b DECLASSIFICATION /DOWNGRADING SCHEDULE 



3 DISTRIBUTION/ AVAILABILITY OF REPORT 

Approved for public release; 
Distribution unlimited 



4. PERFORMING ORGANIZATION REPORT NUMBER(S) 



5 MONITORING ORGANIZATION REPORT NUMBER(S) 



6a NAME OF PERFORMING ORGANIZATION 

Naval Postgraduate School 



6b OFFICE SYMBOL 

e* 






7a NAME OF MONITORING ORGANIZATION 

Naval Postgraduate School 



6c. ADDRESS { City , State , and ZIP Code) 

Monterey, CA 93943-5000 



7b ADDRESS (City, State, and ZIP Code) 

Monterey, CA 93943-5000 



8a NAME OF FUNDING / SPONSORING 
ORGANIZATION 



8b OFFICE SYMBOL 
(If applicable) 



9 PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 



8c. ADDRESS (City, State, and ZIP Code) 



10 SOURCE OF FUNDING NUMBERS 



PROGRAM 


PROJECT 


TASK 


WORK UNIT 


ELEMENT NO 


NO 


NO 


ACCESSION NO 



11 TITLE (Include Security Classification) 

NUMERICAL STUDIES OF LOCALIZED VIBRATING STRUCTURES IN NONLINEAR LATTICES 



12 PERSONAL AUTHOR{S) 



Galvin, Brian R. 



'Wt E e? F ’§ E m 



esis 



13b TIME COVERED 
FROM TO 



14 



DATE OF REPORT (Year, Month, Day) 

March 1991 



15 



PAGE COUNT 

218 



16 SUPPLEMENTARY NOTATION 



The views expressed in this thesis are those of the author and do 
not reflect the official policy or position of the Department of Defense or US government 



17 COSATI CODES 


FIELD 


GROUP 


SUB-GROUP 















18 SUBJECT TERMS ( Continue on reverse if necessary and identify by block number) 

Nonlinear dynamics; Solitons; Nonlinear Stability Theory; 
Chaotic Systems; Numerical Modeling 



19 ABSTRACT (Continue on reverse if necessary and identify by block number) 

A simple numerical model using a modified Euler* s method was developed to model 
nonlinear lattices. This model was used to study the properties of four breather and 
kink type solitons in the cutoff modes of a lattice of linearly coupled oscillators 
with a cubic nonlinearity. These cutoff mode solitons were shown to correspond very 
well to the theoretical predictions of Larraza and Putterman [1984] and the experi- 
mental work of Denardo [1990]. In addition, a fifth soliton was discovered in 
the upper cutoff mode, which was not anticipated by the theory. A preliminary 
analytical attempt to describe this soliton and to describe solitons in the inter- 
mediate modes, due to Larraza, Putterman, and the author, is presented. Additional 
numerical work on intermediate mode solitons and domain walls was performed. These 
studies showed that kink solitons are ubiquitous, and that they appear to be intimately 
linked todomain wall structures. In order to demonstrate the flexibility of the com- 
puter program developed, the model was extended to include two dimensional and 



20 DISTRIBUTION /AVAILABILITY OF ABSTRACT 
^JuNCLASSlF IE D/UN LIMITED □ SAME AS RPT 


□ DTIC USERS 


21 ABSTRACT SECURITY CLASSIFICATION 

UNCLASSIFIED 


22a NAME OF RESPONSIBLE INDIVIDUAL 

Bruce Denardo 


22b TELEPHONE (Include Area Code) 

(408)64§- 


22 c OFFICE SYMBOL 

■ Efl^s 



)D Form 1473, JUN 86 



Previous editions are obsolete 

S/N 0102-LF-014-6603 



SECURITY CLASSIFICATION OF THIS PAGE 



Unclassified 



Unclassif ied 



SECURITY CLASSIFICATION OF THIS PAGE 



lattices and one dimensional lattices with nonuniform characteristics . Two 
dimensional breather and kink solitons are described. Finally, a Toda lattice 
was modeled and some preliminary results obtained in preparation for future 
work. 






DD Form 1473, JUN 86 (Reverse) 



SECURITY CLASSIFICATION OF THIS PAGE 

Unclassif led 



Approved for public release; distribution unlimited 



Numerical Studies of Localized Vibrating Structures 
in Nonlinear Lattices 

by 

Brian Russell Galvin 
Lieutenant, United States Navy 
A.B., Harvard University, 1981 

Submitted in partial fulfillment of the 
requirements for the degree of 

MASTER OF SCIENCE IN ENGINEERING ACOUSTICS 

from the 

NAVAL POSTGRADUATE SCHOOL 
MARCH^.991 



m 



ABSTRACT 



& / 



A simple numerical model using a modified Euler’s method was developed to 
model nonlinear lattices. This model was used to study the properties of four 
breather and kink type solitons in the cutoff modes of a lattice of linearly coupled 
oscillators with a cubic nonlinearity. These cutoff mode solitons were shown to 
correspond very well to the theoretical predictions of Larraza and Putterman [1984] 
and the experimental work of Denardo [1990]. In addition, a fifth soliton was 
discovered in the upper cutoff mode, which was not anticipated by the theory. A 
preliminary analytical attempt to describe this soliton and to describe solitons in the 
intermediate modes, due to Larraza, Putterman, and the author, is presented. 
Additional numerical work on intermediate mode solitons and domain walls was 
performed. These studies showed that kink solitons are ubiquitous, and that they 
appear to be intimately linked to domain wall structures. In order to demonstrate 
the flexibility of the computer program developed, the model was extended to 
include two dimensional lattices and one dimensional lattices with nonuniform 
characteristics. Two dimensional breather and kink solitons are described. Finally, 
a Toda lattice was modeled and some preliminary results obtained in preparation for 
future work. 
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I. INTRODUCTION 



Recent years have seen a great shift in emphasis in physics research as the total 
dominance of quantum mechanical studies has yielded somewhat to the rapidly 
growing group of disciplines commonly grouped under the "nonlinear physics" title. 
What is meant is more specific — the study of classical, nonlinear, dynamical systems. 
Within this group of disciplines is the field of soliton research, which boomed during 
the last three decades. This work has yielded practical results in several areas, most 
notably in the development of fifth generation fiber optic communications systems 
which will transmit soliton pulses, greatly increasing the permissible time-bandwidth 
product of a given system (Hasegawa and Tappert [1973]). Solitons, which are 
spatially localized nonlinear wave packets so named because they have many 
particle-like properties, have became a key part of the current state of the art in 
cosmology, particle physics, condensed matter physics, and hydrodynamics, to name 
but a few. 

Most of the soliton work performed to date has focused on continuous systems, 
such as fluids. Of the little work done on discrete solitons, much has been 
concentrated in the study of cellular automata. Recently, though, a group of 
researchers at UCLA has been conducting research into the characteristics of solitons 
in discrete lattices. Experimental work has been done by Denardo [1990] on a 
simple realization of a nonlinear lattice -- a one dimensional lattice of linearly 
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coupled pendula. He verified the existence of two of four predicted soliton types 
in the cutoff modes of the lattice (the cutoff modes are those modes where all of the 
elements are either in or out of phase with both of their nearest neighbors). 
Additionally, he studied the behavior of these solitons as he varied the frequency 
and amplitude of the parametric drive system he used to overcome the effects of 
damping. 

The purpose of the research presented here was initially to provide numerical 
verification of Denardo’s experimental results, in order to check the degree to which 
his results actually corresponded to theoretical predictions. It should be noted that 
his results were largely qualitative, since the coupling and damping parameters of the 
lattice he used could not be measured directly. A simple computer model of the 
pendulum lattice was developed, tested, and then used to quickly and accurately 
verify the results of Denardo. In fact, the solitons in the pendulum lattice were 
found to be in excellent agreement with the theoretical predictions, which were 
based on a weakly nonlinear approximation. This agreement persisted even outside 
the limits of the approximations used in developing the theory, suggesting that these 
soliton structures are very robust. 

The linear lattice theory basics required for this work are presented in Chapter 
II, along with a discussion of the theory of parametric drive and its effect on the 
stability of lattice dynamics. Chapter III presents the theoretical development upon 
which Denardo’s work was based. Further, the results of a study of the stability of 
various conditions of the cutoff modes is presented. This study was performed in 
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order to determine the effect that discreteness had on the validity of the well known 
stability theory of Benjamin and Feir [1966], which was developed for continuous 
systems. Finally, the numerical results which validated Denardo’s experimental work 
are presented. The existence of the third (which was seen experimentally by Larraza 
et al. [1990] in a water channel experiment), and fourth predicted solitons, and of 
one previously unsuspected soliton were also shown. 

Having quickly achieved the initial goal of validating Denardo’s work, we went 
on to study modes intermediate in frequency and wavelength between the two cutoff 
modes. These results, along with a brief sketch of a theoretical description which 
is still being developed, are presented in Chapter IV. Many different types of 
soliton-like structures were identified and their properties studied in this phase of 
our work, leading us to conclude that, in particular, kink solitons (described in 
Chapter III) appeared to be ubiquitous, existing in every mode we checked. Moving 
beyond simple soliton-like structures, the existence of domain walls was shown for 
many different combinations of modes (these domain walls are essentially localized 
boundaries between two different stable modes). 

Finally, in order to extend the results obtained in the various lattice modes, 
work is presented in Chapter V on the effects of nonuniformities in the lattice on 
the stability and motion of the many different solitons and soliton-like structures. 
Preliminary results from the study of a two-dimensional model of the same type of 
lattice, and from the study of a one dimensional model of the Toda lattice are also 
presented in Chapter V. The Toda lattice is a lattice in which the coupling between 



3 



nearest neighbors is exponential, instead of linear. It has been much studied in the 
two decades since it was first formulated because it has many attractive 
mathematical properties (not the least of which is that it is completely integrable). 

By the time the work presented here was complete, the initial objective had 
been broadened to include not only a study of nonlinear dynamics, but also a 
demonstration of the value of highly interactive computer modeling of real physical 
systems. The numerical method used, which was quite simple, is presented in 
Appendix A; a short manual on using the program and the actual code are presented 
in Appendices B and C. This program, in its final form, was designed to provide a 
fast, accurate model of physical systems which are composed of oscillating elements 
(such as lattices), which was highly interactive and easily modified. The fact that the 
program is highly interactive makes it possible for the researcher to go far beyond 
the traditional method of numerical analysis wherein a set of parameters is put into 
the model and the output is then examined. An analogy can be drawn between the 
traditional numerical analysis methods, which act essentially as filters that alter the 
input in a definite way. The model used in this thesis is analogous to the newer 
adaptive filters, which respond to the dynamics of the system in a way that increases 
the utility of the effort. The adaptive features are provided by the user via keyboard 
commands which are entered as the model is running and which take immediate 
effect. 

The contribution which it is hoped this thesis will make is, therefore, twofold. 
First, many interesting new results were obtained pertaining to solitons in nonlinear 
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lattices. These results have served as a launching point for several interesting 
theoretical developments achieved in our group and at UCLA, and as a means of 
both verifying and directing the work of the experimentalists (whether verification 
or direction occurred depended on who got there first). Second, a simple computer 
model was developed which is highly interactive, somewhat novel in the level of 
freedom it gives the researcher, and easily modified to model other physical systems. 
This program is already being used by several additional students, each of whom is 
making his or her own minor modifications. It is hoped that the program will in the 
end have served as a nucleation site for continuing innovation in interactive physics 
modeling. 
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II. PRELIMINARY: THE LINEAR LATTICE 



A. THE SIMPLE MASS-SPRING LATTICE PROBLEM. 

As a starting point in the study of solitons in nonlinear lattices, we consider 
first the behavior of linearized lattices in this chapter, starting with a simple lattice 
of point masses connected by massless springs. Except in Chapter V, this and all 
other lattices will be taken to consist of identical masses and stiffnesses. Also, 
throughout this thesis periodic boundary conditions are used, such that the last 
element is coupled both to the next to last element and to the first element; this is 
essentially a finite ring lattice, except that all effects of curvature of such a ring 
lattice are ignored. If we restrict our analysis to longitudinal motion, the exact 
equation of motion for the nth element is 

x ~—(x ,+x -lx ), ILA.1 

n ' n+1 n-1 n' 3 

ttl 

with s and m being the stiffness of the springs and mass of the point masses, 
respectively, and with x being the deviation from equilibrium for each element. 
Letting 
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we get 



V“o II.A.3 



We assume an oscillatory solution of the form 




n.A.4 



where na is the spatial variable, and a is the spacing of lattice elements. Substitution 
of this solution into (II.A.3) and solving for frequency yields 



lattice, k is quantized. Thus, as seen in Figure II. 1, the dispersion relation does not 
give a continuous curve such as the one described by II.A.5. 

It is often useful to work with continuum limit approximations of discrete 
systems, so we will derive here the continuum approximation of the linear ring 
lattice. If we let y = na be the spatial variable, and we take the limit of (II.A.3) as 




II.A.5 



which is the dispersion relationship for this linear lattice. For a finite and discrete 
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Figure 11.1. Dispersion relation for a finite ring mass-spring lattice. 
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the spacing of the lattice goes to zero (and the mass and stiffness go to zero and 
infinity, respectively), we get 



a - o a 2 dy 2 



II.A.6 



Thus one can express the equation of motion of the linear ring lattice in the 
continuum limit by a simple wave equation: 

d 2 *" 1 d 2 *" 0 II.A.7 

dy 2 c 2 dt 2 



where 




II.A.8 



Finally, it is worthwhile working backwards for a moment to obtain a result 
that is interesting now and useful later. If we start by letting y be a continuous 
rather than a discrete variable, and consider the value of x(y) at two points spaced 
a small distance a to the left and right of a particular point y 0 , we can express them 
as Taylor series expansions: 
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II.A.9 



-y^K-y. 



and 



M) ,£i 




WU 2 


U 2 J 



+ .... 



n.A.io 



>->0 



Now, adding these together and subtracting twice x(y 0 ) gives us an expression 
identical to the right hand side of (II. A. 1), with the exception of the multiplicative 
constant s, if we ignore the higher order even derivatives. What does this mean 
physically? If we start with a continuous system and discretize it as indicated, we get 
the same equation of motion as our discrete lattice, if we assume only nearest 
nei ghbors interact . Assuming infinite interaction length, for example, means that to 
use (II.A.l) to approximate a continuous system entails an error of order 4 in the 
step size a and in derivatives of x; for harmonic systems, the error is of order (ka) 4 , 
which is small when the second derivative is finite. However, when the second 
derivative vanishes, these errors can be significant. 

B. THE PENDULUM LATTICE AND ITS GENERALIZATION. 

One of the simplest experimental realizations of a nonlinear lattice is a lattice 
of coupled pendula. Such a lattice was used by Denardo [1990] to study solitons in 
discrete lattices. An idealization of the lattice he used is shown in Figure II.2. In 
the actual lattice, coupling was accomplished by tying knots between the V shaped 
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strings which supported adjacent elements, and was assumed to be approximately 
linear. For our purposes, we will assume a completely linear torsional spring 
coupling, for ease of analysis. Each pendulum bob is now acted on by two forces - 
gravity and the coupling force from each nearest neighbor. The exact equation of 
motion is given by 



dt 2 




i - 20„)+o)^sin0 || -O, 



II.B.l 



with 



u> 



o 




II.B.2 



This equation can be rewritten with a coupling constant t replacing the coefficient 
of the second term and the sine term expressed as an infinite series: 



fe. 

dt 2 






U.B.3 



which can then be generalized, with higher order nonlinearities dropped, to 
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Q n -V( e n+l + e n-r 2e n) + “>l e n - ad l 



II.B.4 



In this equation, the sign of a is left unspecified; in Chapter III, it will be seen to be 
of critical importance, whereas the magnitude is not. Throughout this thesis, we will 
take a to be either ± 1 or 0, as desired. 

In order to derive the dispersion relation for (II.B.4), we set a equal to zero, 
since the dispersion relation is a linear concept. As before, we assume a solution of 
the form of (II.A.4), and substitute into (II.B.4). Solving for frequency yields 



G>- 



\ 



o>J+4Ysin 2 |-y 



II.B.5 



This dispersion relation is plotted in Figure II.3. Again, the dispersion relation is a 
discrete function of k, since we are dealing with a finite discrete lattice. As we add 
elements to the lattice, the density of points on the curve becomes greater and 
greater, until the curve is continuous. This corresponds to an infinite lattice. 

Denardo [1990] studied solitons in the two extreme modes represented in 
Figure II.3, referred to as the upper and lower cutoff modes. The lower cutoff mode 
corresponds to k = 0, or to an infinite wavelength, and is the mode where all lattice 
elements oscillate exactly in phase, so that coupling is completely inoperative. The 
frequency of this mode is, as indicated by II.B.5, exactly that of a single linear 
oscillator. The upper cutoff mode corresponds to k = ;r/a, or wavelength equal to 
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Figure 113. Dispersion relation for generalized pendulum lattice. 
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two lattice spacings. In this mode, adjacent elements are exactly out of phase with 
each other. Coupling is maximally effective in this mode, and the frequency is given 
by 



2 2 . 
<*>i-o>o + 4 Y- 



II.B.6 



Taking the continuum limit of the equation of motion exactly as before, we get 



dt 2 



^6 -on 
-v +sm0-O. 

ay 2 



I1.B.7 



Keeping only the leading order nonlinearity and generalizing, we get the nonlinear 
Klein-Gordon equation 



£6 ^ e i ( i ) 2 e _ ae ) H £.8 

dt 2 dy 2 



This is the continuum equation which corresponds to the discrete equation modeled 
in this thesis. As with the mass spring lattice, the model can be viewed as a finite 
element approximation of (II.B.7) with only nearest neighbor interactions allowed, 
or an approximation of (II.B.6) which is accurate, where the second spatial derivative 
does not vanish, to fourth order in ka. 
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C. THE DRIVEN DAMPED LATTICE. 



The experimental work which preceded this thesis was done using a damped 
driven pendulum lattice, and indeed the use of drive and dissipation proved useful 
in the numerical modeling of nonlinear lattices as well. The reason for the need for 
drive in the experimental lattice is that dissipation is unavoidable, so if any steady 
state was to be studied, drive was necessary. In the case of the work presented here, 
the primary utility of drive and dissipation arises from the fact that energy radiated 
away from solitons as they reach a steady state profile is damped and does not 
return after traversing the ring lattice. When the free case is studied, it is often 
difficult to recognize the presence or confirm the absence of solitons in the presence 
of high background radiated energy. 

In the numerical model, drive and dissipation were incorporated into the 
equation of motion by incorporating a linear damping term, with damping constant 
6, and a modification to the last term by the inclusion of a term representing drive 
of amplitude 2r\ and frequency 2o>: 

— Y(*„ +1 + Vi _2 0 + P*„ + [o>J+2Ticos(2tor)]x II - axl UC1 

dt 



This is the equation for a system with acceleration drive , as opposed to a system with 
displacement drive, where we would have an additional factor of o> 2 multiplying the 
cosine in the drive term. Equation (II.C.l), with the coupling term and the cubic 
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nonlinearity neglected, is a form of Mathieu’s Equation, the analysis of which is well 
established (see, for example, Pippard [1979], pp. 289-301, or Appendix I of Denardo 
[1990]). 
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III. THE NONLINEAR LATTICE I: CUTOFF MODE SOLITONS 



A. BRIEF OVERVIEW OF SOLITON THEORY. 

Solitons are "localized regions of motion with finite energy that can appear in 
media that support waves" (Denardo [1990]). They can be viewed either as finite 
extent nonlinear waves or as "particles" (hence their name, which was coined by 
Zabusky and Kruskal [1965]). Solitons have several interesting properties. They 
exist because of the competition between dispersive and nonlinear effects. 
Dispersive effects tend to "spread out" the wave, shedding energy via linear radiation 
of energy at various frequencies. Nonlinear effects tend to concentrate the energy. 
When the medium in question obtains a profile that matches a soliton solution to its 
underlying differential equations of motion, these effects are exactly balanced, and 
energy is trapped in a soliton. Solitons can either be propagating or stationary; 
whereas most of the theoretical and numerical work done until recently has 
concerned itself almost exclusively with propagating solitons, the work here is 
focused almost exclusively on the stationary, or standing, solitons first discerned 
experimentally by Wu et al. [1984] and explained theoretically by Larraza and 
Putterman [1984], Another important aspect of solitons is that they collide elastically 
and are frequently very long-lived. These two facts, and the very great localization 
of energy within a potential field that solitons represent, have made them the subject 
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of much interest in fields as widely separated as molecular biology and particle 
physics. 

Several well established nonlinear equations of motion have been shown to 
have soliton solutions (Ablowitz and Segur [1981]). Of these, we are concerned here 
with only two. The first of these is the Nonlinear Schrodinger (NLS) Equation: 



where, for a surface wave of wavenumber k on a deep liquid, the parameters are 




iu.a.1 




m.A.2 



v-2wjt 2 



m.A.3 



and 




III.A.4 



A well-established soliton solution to the NLS equation is 



A- —sech — (£-vr) e 2y , 



N v N 



III.A.5 



where the velocity v is a free parameter, and where 
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III.A.6 



b-— 

2y 



v ( dm v 



dk 2 



An important feature of this solution is that the amplitude of the envelope divided 
by the characteristic length of the soliton is a constant, with value 



K- 



N 



2y 

v 



III.A.7 



This property proves very useful in discriminating actual soliton solutions and 
solutions which look superficially like the hyperbolic secant soliton but do not obey 
this rule. 

The other standard nonlinear evolution equation with soliton solutions that is 
relevant here is the sine-Gordon equation, 

^- c 2^ +tl) J sin(0 )_o, III. A. 10 

dt 2 dx 2 



which is derived from the equation of motion of the pendulum lattice and in many 
other physical applications (Dodd et al. [1982]). A static kink solution to equation 
(III.A.8) is 
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III. A. 11 



0(x)-4tan \e 



_<-V 

C 



This solution constitutes a transition between two domains separated by 360 degrees, 
similar to the separatrix motion of a single pendulum. Solutions to the nonlinear 
evolution equations (III.A.l and III.A.8) have been shown to strictly meet the 
accepted standards of what a soliton is. That is, they collide elastically and have 
finite spatial extent. As the exact equations of motion of the generalized nonlinear 
lattice are explored in subsequent sections, it will become clear that these are only 
approximate solutions to the real lattice. No attempt has been made in this thesis 
to examine the elasticity of collisions between two or more solitons, so the term is 
used somewhat more loosely here than in the literature. The width-amplitude ratio 
test has been used, however, to limit the study to solutions at least very closely 
resembling solitons. In fact, it is a general observation that, in any real systems, 
which therefore have some dissipative effects, however small, none of the exact 
soliton solutions that have been discussed in the literature to date actually remain. 
All of those studies have been in free systems, since a dissipative term in the 
nonlinear evolution equations invariably separates the equation form the standard 
forms which have exact soliton solutions. 
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B. NLS THEORY OF THE CUTOFF MODE LATTICE. 

We consider in this thesis the generalized equation resulting from a linearly 
coupled lattice of nonlinear oscillators. The original system in question was the 
model pendulum lattice whose exact equation of motion is 

11181 

where 




III.B.2 



and 



Y 




III.B.3 



and T is the torsional constant of the springs connecting adjacent elements, a is the 
lattice spacing (which we henceforth take to be unity, and drop from the equation), 
m is the mass of the pendulums bobs, and L is the length of each pendulum. This 
equation is usually linearized as 

0 -v(0 ,+0 ,-26 )-Wn6 . III.B.4 

n /i+l /*- 1 n' w O v jj- 

However, in the weakly nonlinear regime, the second term in the Taylor expansion 
of the sine function is included: 
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III.B.5 



(0 



®.-Y(e„. 1 *e,.,-2e 11 )*»je lt ~^e 3 , 



Generalizing this equation as was done in Chapter II, and adding drive and 
dissipation, we get the equation used in the computer model used in this thesis (for 
description of numerical methods used, see Appendix A): 



0»-Y( e „.r e ..r 2e .)-P®»-(“o’- 2 '> cos ( 2 “')) 8 »- ,ae «' In - B - 6 



Analytically, this equation is difficult to deal with. We consider first, therefore, 
only the two cutoff modes (see Chapter II), with linear frequencies <i> 0 and o>„ 
respectively, where it is possible to make suitable approximations and arrive at an 
analytically tractable expression (in fact, the NLS equation). We consider the 
uniform lower cutoff state to be modulated by an envelope that is slowly varying in 
both space and time. That is, 

e(x,0-A(x,0e' o, +B (x^y 3 ^‘+...+c.c., III.B.7 



where 

tBWAi III.B.8 

and the lack of even harmonics is due to the cubic nonlinearity. Substituting this 
solution into the cubic Klein-Gordon (III.B.6), dropping all higher harmonics, and 
neglecting A tt compared to A, because A(x,t) is slowly varying, we get 
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dt dx 



III.B.9 



for the lower cutoff mode and 



2/o>— +Y-^^+(o)i-o> 2 +/o>P)A+riA-3alAl 2 A, 



dt dx 2 
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for the upper cutoff mode, recalling that 




m.B.ii 



Thus both the upper and lower cutoff modes can be represented by the nonlinear 
Schrodinger equation, in the continuum limit, and at the weakly nonlinear level. The 
next theoretical task is to examine the stability of the uniform states of these modes; 



including their stability characteristics. 

C. STABILITY OF THE CUTOFF MODES. 

Following the method of Stuart and Diprima [1978] as extended by Denardo 
[1990], we use an amplitude equation to study the stability of the uniform cutoff 
modes of our generalized lattice. Stewart and Diprima showed that this method is 
equivalent to the somewhat more complicated perturbative methods used by 
Benjamin and F ir [1967] and Eckhaus [1965]. The basic idea is to consider a small 
perturbation in the sidebands of the uniform mode and to determine whether the 



following that, the soliton solutions of the NLS equation will be examined in detail, 
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mode is stable when so perturbed. The essential result of Benjamin and Feir, and 
of Eckhaus in a slightly different but conceptually equivalent problem, was that, for 
certain wavelength regimes, resonances between the sidebands and the mode caused 
energy to be shed by the mode to the sidebands, which grow exponentially. 

Eckhaus and Benjamin and Feir were concerned with fluid phenomena when 
they conducted their work. In fact, their results were developed entirely in the realm 
of continuum mechanics. However, they can be used profitably with continuum limit 
approximations of discrete lattices, such as the Nonlinear Schrodinger Equation. 
However, the effects of finite lattice size on the theoretical thresholds of instability 
have not been studied heretofore. Later in this section we will consider the case of 
a two element oscillator lattice under the equations of motion given above. It will 
be shown that one can solve approximately for the threshold of instability in this 
limiting case, and that the result is in fact different from that predicted by the 
continuum theory (which is not really surprising; rather, it would have been 
surprising if the results had been identical). Finally, numerical results for lattices 
from two to 40 elements in size will give a very clear relationship between the 
continuum theory and the actual finite lattice results. There currently exists no 
theoretical framework in which to place these results, other than to compare them 
with the two element lattice theory (with which the agreement is remarkable). 

1. Free Lattice Stability. 

We consider the NLS equation for the lower cutoff mode of the free 

lattice: 
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2;o)i^.-c 2 — + (o)J-(o 2 \A-3alAFA. 
dt ax 2 [ r 



The uniform steady state is A = a positive real constant, with the relationship 
between amplitude and frequency being 




m.c .2 



To investigate the stability of the motion, we set 



A“A 0 (1+T), 
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where 



Y=Y(x,f) 



III.C.4 



and 



lYkl. 



m.c.5 



Substituting this into the NLS equation, neglecting higher order terms in Y, and using 
(III.C.3) to eliminate zero order terms, we get the evolution equation (Denardo 
[1990]) 
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2;wT r c 2 T i;r -3a/lo(T + T)-0. III.C.6 

Imposing a small modulation of wavenumber k, and including both left and right 
traveling waves, we set 

y-ae^'Kbe*-*^, llLCJ 



with a, b, and Q complex constants. Grouping terms from (III.C.26) according to 
the type of exponential factor, and setting the coefficient of each to zero, we get 



2o>Q+c 2 fc 2 -3aj4o 


-3aAl 


a 


-3aA 2 0 


-2(*>Cl+c 2 k 2 -3aAo 


b 



III.C.8 



Setting the determinant of (III.C.28) to zero and solving for Q (Denardo [1990]) 
gives 



Q- 



c 2 k 2 



2(o 



1- 



6olAq 

T 2 * 2 " 



III.C.9 



Thus, for Q to be real and the modulation to be stable, we require 
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6aAo<c 2 k 2 . 



Thus, for the lower cutoff mode, if a>0 (softening lattice), the motion will be 
unstable at a given amplitude to a sufficiently long modulation wavelength. Since 
the wavelength is limited by lattice size for finite lattices, it is clear that for finite 
softening lattices there exist amplitude thresholds below which the lower cutoff mode 
is stable under modulation. For hardening lattices, an identical procedure starting 
from (IU.B.10) yields 

6a A Q>c 2 k 2 . HI-C.11 

Thus, for a hardening lattice, the free uniform state will be unstable for a given 
amplitude for modulation wavelengths beyond some threshold, exactly analogously 
to the softening lower cutoff mode. Also, the softening upper cutoff mode is always 
stable to sideband modulations. This result is in fact the first of many which show 
the marked symmetry properties of the cutoff modes. 

2. Stability of the Uniform Damped, Driven Cutoff Lattice. 

Starting from (III.B.9), and following the method again of Denardo [1990], 
we consider first the effects of parametric excitation on the nonlinear cutoff modes, 
then the stability of the damped driven lattice to sideband modulation. Noting that 
(in.B.9) is equivalent to (II.C.4), with the addition only of the cubic nonlinearity 
term, it is clear that the analysis of excitation from rest is unchanged from that 
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performed in Chapter II, since, at and just above the state of rest, the additional 
nonlinear term plays no role. However, as mentioned in that analysis, the 
nonlinearity does, above some amplitude, begin to play a role that counteracts the 
exponential growth of oscillation due to the parametric excitation. Therefore, where 
the linear case grows without bound once the excitation from rest threshold is 
reached, the nonlinear lattice does not. 

Consider the stability of these damped driven lattice states when 
modulated by sidebands. We start as before with 



To investigate continuum stability for a modulation of wavenumber k, we let 



As before, we retain only those terms linear in Y, which gives us (Denardo [1990]) 




dt dx 



A-V^(1+TW), 
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where 
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III.C.14 
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and 
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Letting 



y-ae jikx *' 1 ' ) +be' i(kx *' 1, \ 



III.C.19 



we get the coefficient matrix equation 

f-htf-v. _ g |«1 m.C.20 



We thus find that the motion is stable if 
and 
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i-hk f*P, 
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which, when the values for f and h are substituted in, yields 



3aAo~-^c 2 k 2 



! |^3aylo-^ 2 /: 2 j-(coo-o 2 ) >0. 



III.C.23 



This result, from Denardo [1990], can be expressed graphically for an acceleration 
drive system, as shown in Figure III. 1. To use Figure III.l, select any point on the 
tuning curve. Decrease the ordinate by (^1^/2. If the resultant point lies in a cross- 
hatched region, then the state is unstable to a modulation of wavenumber k. 
Similarly, Figure III.2 shows the stability region for the upper cutoff mode, which is 
essentially the same as the lower cutoff mode, with coj replacing ti> 0 . 

3. Effects of Finite Lattice Size on Stability. 



continuum theory based on the work of Benjamin and Feir [1967]. However, in 
experimental and numerical work, finite lattice sizes are a necessity, even when 
periodic boundary conditions are used, as they were in all of the work presented 
here. Theoretically, although it is clear that finite lattice size leads to a discretization 
of possible modulation wavenumbers, the effects of finite lattice size on stability 
remain incompletely understood. In an effort to remedy this situation, we consider 



The results presented in sections 1 and 2 were based entirely on a 
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the case of a two element lattice, which can be solved more directly (although still 
only approximately) using the equations of motion. This analysis will be performed 
and the results compared directly to those of the previous Benjamin-Feir analysis. 
Then numerical work on the subject will be presented, which connects in a 
continuous fashion the results of the Benjamin-Feir analysis and our straightforward 
N = 2 analysis. 



a. 



Two Oscillator Lattice Stability. 



For the two oscillator lattice, the exact equations of motion are 



e 1 +wj0 1 -ae5--Y(e 1 -0 2 ) 
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e 2 +G)oe 2 -a02-Y(e,-0 2 ). 
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Transforming to normal coordinates according to the transformation 






III.C.26 



C-i(8,-9 2 ), 



UI.C.27 



we get, after substitution into the equations of motion, 



and 
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Some manipulation shows that 



e,±e 3 - 



2Z(Z 2 +3( 2 ) 

2C(C 2 +31; 2 ) 
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This then yields, for equations of motion in normal coordinates, 



l+ult-aZ 3 + 3aC 2 Z 
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and 



C+(o)o + 2y)C-aC 3+ 3aS 2 C- 
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An examination of the form of these equation shows that each of the 
modes, which can be considered as oscillators, is parametrically driving the other. 
We need to determine the amplitude threshold for excitation next. Consider, for 
weakly nonlinear motion, the initial state 



C -Acos<i3t+(higher harmonics), 



III.C.33 
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where 




in.C.35 



For a single parametrically driven oscillator, 



0 +(Q 2 +2ti cos2g)00 -0, 



III.C.36 



the amplitude threshold is given approximately by 



t|>|Q 2 -o) 2 |. III.C.37 

This is only valid if the two frequencies are very similar, and is often expressed as 

r|>2QlQ-6)L III.C.38 



We neglect the cubic term in (III.C.31) because the amplitude is small initially. The 



other nonlinear term is 



3aA 2 cos 2 o>t^ - — A 2 ( 1 +cos2o>r) £ . 
2 



III.C.39 



There is, in addition to the driving force, a restoring or antirestoring force which 
alters the frequency of the driven oscillator. The new frequency is given by 



Q 2 -g>o~— Oj 4 2 . 
0 2 
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The amplitude of the parametric drive is given by 

2rj — — laU 2 . 

2 
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Excitation occurs if 



— Ial4 2 >| Wq- — co4 2 j - ( g)q + 2 y — co4 2 
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The alteration of the frequency of the driven oscillator is very important since it is 
a greater change than the nonlinear bending of the frequency of the driver oscillator: 
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This relation is not satisfied if a > 0 (softening case), so the softening upper cutoff 
lattice is stable (as expected from previous considerations). For a < 0, it follows 
from (III.C.42) that the threshold condition is 

— laL4 2 >y. III.C.44 



Therefore we can define the critical amplitude to be 



N 



4y 

3lal 
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A similar analysis shows that, for the lower cutoff mode, the motion is unstable for 
amplitudes greater than Before we can compare this critical amplitude, which 
applies only for the N = 2 case, to the Benjamin-Feir critical amplitude (obtained 
from III.C.10), 
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we need to express it in terms of the original coordinates (Equation (III.C.45) is for 
the normal coordinate system): 



A 



C//-2 
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For a coupling parameter of 0.25, the values are 1.28255 and 0.816497, respectively. 
It is clear that the difference is significant between the continuum theory prediction 
and the discrete theory; any attempt to explain theoretically the deviation from 
continuum theory in the finite lattice limit must account for this difference. 

b. Numerical Study of Finite Lattice Effects. 

A simple lattice model, described in Appendix A, was used to explore 
the behavior of the stability threshold for finite lattices in order to explore more 
deeply the relationship between Benjamin-Feir continuum theory and the exact 
theory just presented. In each case, a lattice of amplitude A was modulated by a 
small (about 10% peak to peak of A) modulation of wavelength equal to lattice 
length. The long term behavior of the system was then monitored. If the initial 
amplitude A was above the actual stability threshold, then the clean, single spatial 
frequency modulation grew and became complex spectrally. Otherwise, the system 
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remained stable for long (at least 100 cycles). Various amplitudes were used, with 
increasingly fine resolution, until a good estimate of A,, was available. All values 
shown in graphs are accurate to within less than one percent. 

First we consider the N = 2 lattice as a special, limiting case of the 
finite lattices. Figure III.3 shows the theoretical curves of A,, as a function of 
coupling, along with the results obtained for a hardening lattice with numerical data 
obtained using the model. It is clear from the numerical results that the approximate 
theory is correct, which was expected. Extension of the exact two element theory 
to larger lattices has not yet been attempted, but a numerical determination of the 
onset of instability as a function of lattice size was conducted. The results are given 
in Figure III.4. The continuum theory clearly determines the threshold for lattices 
larger than 20 elements; for smaller lattices, a mechanism probably related to the 
one studied in the two element case causes instability to occur sooner. 

While verifying this theory and showing the extent of deviation from 
continuum theory that the smallest possible finite lattice exhibits, an interesting and 
unexpected additional result was obtained. Before describing it, it is important to 
point out that the instability predicted by the exact theory in the last section was due 
to a different mechanism than the Benjamin-Feir instability. Instead of sideband 
modulational instability, the exact theory merely applied well established (and 
previously discussed) parametric excitation stability theory to the case where one 
element is considered to be parametrically driven by the second element. In fact, 
one might have been tempted to presuppose that the results obtained using the 
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Figure III.3. Comparison of theory and numerical results for N = 2 stability. 
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Fig. 111.4. Numerical stability results, y = 0.25, with Benjamin-Feir theory. 
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parametric drive theory would be invalid, since, as the driven element’s amplitude 
grows, it begins to be stabilized by its nonlinearity. Additionally, at some point one 
would guess that, in a phenomenon similar in appearance to the well understood 
linear energy beating phenomenon, the driven element becomes the driving element. 

In fact, the instability which occurs at the threshold amplitudes 
calculated by the exact parametric theory is a weak instability . For amplitudes just 
slightly greater than the threshold, the driving element periodically decreases and 
then increases in amplitude, while the driven element does the opposite. For 
amplitudes in the lower third of the band between the two theoretical thresholds, the 
two elements never switch roles. Then, there is a band where the roles switch, but 
in a chaotic fashion. That is, the switching takes place at unpredictable intervals. 
At slightly greater amplitudes, the switching becomes more periodic, but the detailed 
motion of each element becomes more and more chaotic. The boundary between 
the domains where switching does and doesn’t take place proves to be very difficult 
to pin down, as in fact do most of the similar boundaries encountered in this study. 
This can be seen, for example, in the time series of Figures IU.5 and 111.6. Each of 
these is a series of measurements of the amplitude of the first of two oscillators in 
a lattice with coupling of 0.15. In Figure 5, the lattice started with amplitude of 
0.8225; in Figure 6, the initial amplitude was 0.84. For reference, the lower 
threshold amplitude for this coupling is 0.65, and the upper threshold is 0.99. In 
Figure 5, the lattice is just inside the boundary, that is, it is just inside the no 
switching domain basin of attraction. The lattice makes a sudden flip across the 
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boundary, and then flips back again. This behavior occurred in a free system after 
the system had been running for several hundred oscillator cycles. In Figure 6, the 
lattice lies initially just outside the boundary, then flips back and forth a couple of 
times before settling inside the boundary where it stayed for several thousand 
periods of the oscillator. The entire transient took dozens of oscillator cycles, but 
only fifteen of the long time cycles of the drive. Based on Figures 5 and 6, one can 
conclude that the boundary lies somewhere around 0.83. A more striking example 
of a state lying exactly on the boundary is given in Figure 111.7. The element is from 
an N=2 lattice with maximum (0.25) coupling and initial amplitude of 1.05, which 
is 51% into the inter-threshold band. The extreme irregularity and clear lack of a 
preferred state of this free system are typical of highly chaotic systems. 

While watching the motion in the time domain gives a vivid example 
of chaotic motion, it is in the frequency and phase domains that the best careful look 
at the onset of chaos can be had. In order to fully characterize the evolution of the 
system’s behavior between the two theoretical threshold amplitudes, a study of 
behavior at various amplitudes was conducted using frequency and phase domain 
methods . Figures IU.8 and II1.9 give a sampling of the spectra of one of the 
elements at various amplitudes. It is important to note that these are representative 
samples only, for the spectra for most amplitudes are continually shifting, especially 
in the higher amplitudes where there is frequent role shifting between driver and 
driven elements. 
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Figure III.7. Chaotic N=2 lattice element time series with max 




Figure 111.7 (cont.) 
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Examining the spectra in Figures III.8 and III.9, one sees that the 
element is monotonic just below the lower threshold (the exact theory threshold). 
As one moves to just above the threshold, several small tonals appear very close to 
the base frequency, and the signal becomes somewhat more broadband in character. 
At 20% into the "inter-threshold band", the side tonals have become more 
pronounced, and have moved farther from the home frequency. Note also that they 
are mainly on the high side of the main tonal -- for amplitudes in the weakly 
unstable band, this is typical of the driven elements. The driving element typically 
has sidebands on the low side of the main tonal. At 51% into the band, a great deal 
of broadband noise is seen; motion in this regime is chaotic. The number of 
sidebands is also much greater. Notice, in the first example (peak amplitude 652) 
that most of the tonals are on the low side, and that the spectrum is somewhat 
cleaner than the other. This is indeed a driving element, although the roles switch 
in this portion of the band. The second spectrum (peak amplitude 982) is typical of 
an element that is beginning to switch. A relatively clean spectrum becomes very 
noisy, or chaotic, and then it switches from predominantly low to high tonals, or vice 
versa. This corresponds to the kind of chaotic transient seen in Figures 5, 6 and 7 
in the time domain, and is in fact for the same conditions as those which gave the 
markedly chaotic time series in Figure IH.7. 

The next pair of spectra are at 72% into the inter-threshold band. 
The first is remarkably clean, and corresponds to a driven element, even though 
there is an imbalance to the low side of the main tonal. In fact, as one moves higher 
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Figure 111.8. Spectra, N = 2 lattice elements just below and above lower threshold. 
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Figure III.9. Spectra for N=2 lattice elements at 51% and 75% of unstable band. 
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into the band, one sees the distinction in sideband distribution between driver and 
driven element break down, and one has to rely primarily on peak amplitude to 
classify. The second element, based on this, is clearly a driving element, and it still 
exhibits the predominant low end sidebands. It is also beginning to switch, as 
indicated by the greater amount of noise. 

The trends seen in this small sampling of spectra in the inter- 
threshold band are valid throughout the band. As one moves from just below the 
lower threshold to just below the upper threshold, one sees a gradual move from 
monofrequency response to a response with many sidebands. The number of 
sidebands grows as amplitude increases, but it does not grow with time at a given 
amplitude (this is different than what one would expect with a Benjamin-Feir type 
instability). There is a clear distinction low in the band between driver and driven 
element, both in amplitude and sideband location. This distinction due to sideband 
location breaks down as one moves above the midpoint of the band, due to the 
element spending an equal amount of time in each of the roles. As seen above in 
Figures III.5 and III.6, an element near the boundary spends large segments of time 
in one role or the other, then may switch in a chaotic burst of activity to the other 
mode. Regular switching of roles between driver and driven does not begin until 
amplitude is well above the boundary between the two domains of attraction. At 
first irregular, or chaotic, it becomes regular as one moves up in the band, which 
corresponds to what has already been seen in the time domain. 
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Finally, looking at the lattice between the two amplitude thresholds 
in the phase domain yields additional insight into the behavior of this complex set 
of states. Figures III. 10 and III. 11 show a succession of phase diagrams for the same 
lattice element shown in Figures 111.8 and III.9 in the frequency domain, and Figure 
III.7 in the time domain. The values of amplitude used are also identical to those 
of Figures III.8 and UI.9, allowing direct comparison. The first phase diagram, 
corresponding to an amplitude just below the lower threshold, is very clean -- this 
is stable, monochromatic motion, as seen in the spectrum. In the next diagram, 
taken for amplitude slightly above the threshold, the same underlying shape is clearly 
evident and takes the majority of the "hits" of the Poincare section, but there is now 
a significant amount of scatter both in and out of the original shape. The idea of 
"weak instability" is here given a vivid graphic representation. 

The phase diagram for 51%, which corresponds to the highly chaotic 
condition at the boundary discussed above, the original shape is still visible, but its 
banded structure has broken down. The region is more amorphous, although some 
hint of the banded structure can still be discerned. The amount of scatter is greatly 
increased, and there is a slight darkening in the inner part of the ellipse, 
corresponding to the fact that, with switching of roles taking place, if on an irregular 
basis, does lead the element to spend an above average amount of time at small 
energy levels. In fact, the well defined outer boundary of the ellipse corresponds to 
the maximum energy state, that is, the state wherein all of the system’s energy is 
concentrated for a brief moment in this particular element. Points outside this ellipse 
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Figure III.ll. Phase diagrams for 51% and 75% of band. 
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are not possible in a free system with the initial conditions given. The fact that this 
boundary is seen here and not in the previous phase diagram is due to the failure, 
low in the weakly unstable band, of the driving element to drive the other until it 
is completely spent (at which point the other element has all or nearly all of the 
system energy). Switching of roles cannot take place until this outer ellipse is 
reached at least occasionally; with the ellipse so well defined, the existence of 
switching, even chaotic, is reasonable. At 75% of band, the second diagram in 
Figure III. 11, the regular switching of the system is readily apparent from the 
existence of the dark orb in the center and the dark ellipse without structure that 
exists halfway out to the maximum energy ellipse. Thus the use of phase diagrams 
corroborates and illuminates what has been learned from the time and frequency 
domain results. 

Figures III. 12, III. 13 and III. 14 show the behavior of lattice elements 
just above the Benjamin-Feir threshold. In Figure III. 12, the time series of the first 
element of a lattice with 0.15 coupling is seen to exhibit regular switching, but the 
period of the switching is variable about an average value. When compared to the 
series in Figure IH.5, which used an identical time step, it can be seen that the 
average period is roughly the same; in fact, this value appears to be a function of 
coupling only. One could visualize the envelope of this regular switching behavior 
as simply a dnoidal wave (Figure m.5) with kinks between every peak, although this 
is rather a bold step of imagination. Figure III. 13 shows various spectra for the 
maximum coupling case (the same lattice as in Figures IU.3 and m.4). The third 
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Figure 111.12. Time series, N = 2 lattice just above Benjamin-Feir threshold. 
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spectrum was the last taken; the growth of sidebands is clearly evident. Finally, 
Figure III. 14 shows the very diffuse phase diagram. One can see readily, as in 
Figure III. 11, that switching is occurring, but there is little else to glean from the 
phase diagram. 

This rather detailed look at the band between the approximate N = 2 
theory threshold of instability and the Benjamin-Feir continuum theory of instability 
support the idea that the two theories are complementary. When the first threshold 
is reached, a weak instability begins to manifest itself due to the passing of the 
parametric excitation threshold amplitude. As amplitude increases, a chaotic 
boundary region is traversed wherein switching between driving and driven states 
occurs irregularly; this boundary region passes into a region where switching is 
regularly, although not precisely periodic. Above this, one reaches the Benjamin- 
Feir threshold, where the mechanism for strong instability is the dumping of energy 
into successive sidebands. The transition between the two types of instability is not 
clear or dramatic, but it is apparently real. 
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Figure III. 13. Spectra for element just above Benjamin-Feir threshold. 
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Fig. III. 14. Phase diagram for an element just above Benjamin-Feir threshold. 
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D. THEORY OF CUTOFF MODE SOLITONS. 



1. NLS Solitons in the Cutoff Lattice. 

Larraza and Putterman [1984] showed that, for the static case, with 
frequency below the linear lower cutoff frequency and a softening lattice, a breather 
soliton is a solution: 
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When the sign of the nonlinearity is changed to hardening, a solution was shown to 
be, in the p = 1 limit: 
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(Denardo [1990]). When considering the upper cutoff case, Denardo showed that 
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is a kink soliton solution of the NLS equation for a softening lattice, if ti) 0 is replaced 
by o), and 
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is a breather for the hardening lattice. 

With the exception of this last solution, Denardo [1990] observed all of the 
soliton solutions for the NLS approximation of the discrete lattice experimentally. 
However, for the pendulum model results, the observations were only qualitative, 
since it was not possible to accurately measure the intrinsic lattice parameters of 
damping and coupling, and the coupling itself was probably not truly linear. This in 
fact was the starting point that motivated the work of this thesis -- to verify the 
experimental results numerically and compare them with the theoretical predictions. 
This work is presented in the next section. 

Before we proceed to the numerical work, some observations on the 
symmetry properties of the theoretical predictions and their underlying physical basis 
need to be made. The symmetry is dual for the cutoff modes -- there is symmetry 
about zero in a, and there is symmetry from upper to lower cutoff. We note by 
examining the equations for the various soliton solutions listed that the parameters 
are identical for both breathers, with the exception of substitution of frequencies and 
sign changes which are necessary because of the changes in a and in position on the 
dispersion curve. The same holds true for the kink solutions. 
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E. NUMERICAL OBSERVATIONS OF SOLITONS IN CUTOFF MODES. 



In his experimental work, Denardo [1990] showed that three of the four 
solitons predicted by the NLS theory for cutoff lattices did indeed exist. He was 
unable to confirm the existence of one, the upper cutoff hardening breather. In 
addition, he mapped the drive planes of the solitons he observed. These mappings 
consisted of varying drive frequency and amplitude and determining the region in 
the plane so mapped in which a given soliton was stable. His lattice work 
necessarily suffered from a lack of quantitative comparison to theory, however, 
inasmuch as the actual coupling and damping parameters of the lattice were not 
determined. This was the initial starting point of the research presented here -- to 
provide quantitative corroboration of Denardo’s results and to compare them with 
the NLS theory. This effort in turn led to many additional areas of exploration, 
which are detailed in later chapters. 

Numerically, all four of the NLS solitons were found to exist in stable states. 
An additional soliton-like structure, a hardening upper cutoff kink, was also 
identified. In the four NLS cases, the agreement to theory was found to be 
excellent, with the degree of departure from theory apparently diminishing linearly 
as the time step is decreased. At a time step of one percent of a period, the error 
was less than one percent, and it got better from there (see Appendix A). This close 
match to theory also naturally led to the width-amplitude product being constant for 
given lattice parameters, as indeed it was. 
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In the detailed discussion of lattice model results in this and later chapters, we 
consider scaled displacements such that the magnitude of the nonlinear coefficient 
is unity, with only the sign varying from hardening to softening systems. The natural 
frequency for a single uncoupled oscillator is also always set to unity, except in 
Chapter V, where nonuniformities are discussed. Thus, the lattice in each case is 
characterized by coupling and damping constants. The drive is characterized by its 
amplitude and frequency. Thus, the parameter space is four dimensional; however, 
only a small part was investigated here. It is also important to note early in the 
discussion of actual lattice results that they are not unique, for a given point in four 
dimensional parameter space. For a given point, there can exist no solitons, there 
may exist one or more soliton solutions, and the solitons that do exist may exist in 
groups or singly. Thus initial conditions play a critical role in the ability or inability 
to achieve a desired state. Also, since the program used was highly interactive, it 
was possible to approach a point in parameter space in several different ways, each 
having slightly different results. As an examination of the tuning curves for 
parametrically driven oscillators makes clear, for example, the system is much more 
sensitive in general to drive frequency changes than to drive amplitude changes. 

To give a detailed example of the challenge offered by the complexity of the 
problem, consider the attempt to reach a state that is in the upper left corner of the 
drive plane region of stability for a structure -- that is, a state where drive frequency 
is low and drive amplitude high. If one tries to reach this state from a precursor 
state that is high in amplitude but at a higher frequency, the transition will be 
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impractical or impossible (I do not know which), since, for any reasonable frequency 
increment, the transient at high drive amplitude is far too great, and the system 
invariably either blows up or transforms into an unrelated lattice structure. If, on 
the other hand, one approached the state by starting from a higher frequency, lower 
amplitude point and moving frequency to the desired value the transients will be 
manageable, and one can then easily increase the drive amplitude to achieve the 
desired state. This simple example points out the characteristics of this work which 
made it very much like trying to chart out an unexplored sea. It also emphasized the 
great value of interactive numerical analysis. Very few of the results achieved in this 
thesis would have been possible if a traditional numerical routine where one puts in 
a parameter choice and gets an output, then repeats the process, is used. 

We begin our look at cutoff solitons with the NLS cases, shown in Figures 
III. 15 through III. 18. The upper and lower cutoff solitons are, theoretically, 
symmetrical, as discussed before, so we will analyze in detail only the hardening 
cases. They offer a richer harvest, since they are not constrained by a potential 
energy barrier, and so can exist at arbitrarily high amplitudes. The theoretical plot 
for each is shown, and the good agreement is apparent. The evident difference at 
the low amplitude elements in the kinks are characteristic of all of the numerics -- 
for a given time step, the overshoot error is larger for smaller amplitudes. As in the 
larger amplitudes, however, the smaller amplitude elements can be brought 
arbitrarily close to theory by shrinking the time increment sufficiently. 
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*igure III.15. Lower cutoff hardening kink example. 
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Figure III.16. Lower cutoff softening breather example. 
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Figure III. 17. Softening upper cutoff kink example. 
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Figure III.18. Hardening upper cutoff kink example. 
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While these solitons were somewhat difficult to achieve initially, they are 
remarkably robust structures. In fact, in the case of the kink, the only way to get rid 
of it is to annihilate it with an antikink (which was done numerically on numerous 
occasions), or to destroy the cutoff mode. When one increases or decreases the drive 
amplitude, the kinks sharpen or broaden accordingly, while maintaining there 
position and identity. However, there is one way that one can force a significant 
change in the kinks. As the amplitude is increased, the kink sharpens until it gets 
to a point where the structure with a node element is unstable to random 
perturbations. The kink undergoes a transition, and arrives at a state where the 
node is between two elements, as in Figure III. 19. This is effectively a one half 
lattice site shift in the position of the kink. This phenomenon, which will be seen 
again in Chapter IV, is not well understood theoretically; in fact, as the studies in 
Chapter IV will make clear, there is in some circumstances a definite preference in 
the discrete lattice for the masses to be at certain definite points of the underlying 
waveform (i.e., the nodes or the peaks, etc.). An interesting phenomenon of the 
breather, noted experimentally by Denardo, is the existence in the drive plane 
(Figure III.20) of a region where a quasiperiodic state exists (Denardo [1990]). 
Further, if one proceeds through this quasiperiodic region, a region where a highly 
self -focused state exists is reached. This state was characterized in the experimental 
lattice by one element oscillating with an amplitude of about 50 degrees while the 
two adjacent moved slightly and the remainder of the lattice was at rest. This 
behavior is remarkable for a system that is driven by a global parametric drive - 
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Figure III.19. Hardening lower cutoff kink with offset node. 
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attempt was made to map the drive plane to the extent done experimentally by 
Denardo (work which will be done by a follow on student), the existence of these 
two states was verified quite early on numerically. Figure III.21 shows the self- 
focused hardening breather. Since this lattice is numerical and not subject to the 
realities of physical systems, the ratio of amplitudes between the main and adjacent 
elements is remarkably greater than 1000! 

The quasiperiodic state is interesting in its own right, of course. The period is 
observed to be stable in these systems; the amplitude of the large element grows and 
decays in a regular cycle. These observations correspond to experimental 
observations made earlier by Denardo [1990]. This quasiperiodic state can be 
understood, I believe, as a transition zone similar in general concept (although not 
in detail) to the transition zone between kinks with and without nodes at masses 
(there is a range of parameters where either state can exist; it is not until one leaves 
this range that the transition takes place from one to the other). To the right and 
below the quasiperiodic region in the drive plane , the breather essentially occupies 
five lattice sites. Above and to the left of the quasiperiodic region, it occupies only 
three sites. In the quasiperiodic region, it seems that the two elements adjacent to 
the high amplitude element first drive the fourth and fifth elements, and are in turn 
driven by the main element. At some point, they cease to drive the outer elements, 
which then drive the inner elements until they decay and the process repeats. Thus 
the quasiperiodicity may be viewed as a nonlinear beating process similar in 
character to the case studied in the two element lattice in Section 3. These self- 



72 



4 . 2 



3.0 





Figure III.20. Drive plane for upper cutoff breather (from Denardo [1990]) 
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links, and may in some similar form be significant as biological solitons, although 
these connections are speculative. 

The fifth type of cutoff soliton was discovered quite by accident. It is a 
"positive energy" hardening upper cutoff kink, shown in Figure HI.22. By positive 
energy, we refer to the fact that the total energy of a lattice with this kink is greater 
than that of a uniform lattice (the NLS kinks are by contrast negative energy kinks; 
in fact, the ability exists with the lattice program to measure the energy of kinks and 
other structures. This was not done due to time constraints, and may be an item of 
interest in follow on work). Due to the late date of their discovery, only preliminary 
study of them has been conducted. They are effectively "brightons", which are 
solutions of the form 



A(x) -K s jUasech\Z(x-x ( )) UI.E.l 

where K and Z are dependent on lattice parameters (Larraza and Putterman [1991]). 
Exact expressions for K and Z have not yet been derived; this theory is for 
conceptual purposes. Additionally, these brightons appear to involve violations of 
the Benjamin-Feir stability criteria; this will be explored in more detail later. 

The brightons, like the rest of the cutoff mode solitons, are remarkably robust 
structures. They maintained their identity even when the time step was increased 
to 25 percent of a period and a perturbation of five percent was added. 
Topologically, the smallest brighton is one in which exactly two elements are 
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Figure 111.21. Hardening upper cutoff self-focused state. 
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! "igure III.22. Low amplitude hardening upper cutoff brighton. 
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positioning occurred at lower amplitudes . As the drive amplitude was lowered to 
0.1, the structure of Figure III.22 broke up and underwent a prolonged and confused 
transient, ending up in the state shown in Figure III.23. This is a very complex state. 
The brighton appears to have become a darkon, but is equally clearly not an NLS 
hyperbolic tangent kink. The mode itself has undergone a transition into something 
which preserves wavelength at two but is otherwise quite different. It is not merely 
shifted spatially, for no shift could give the amplitude division shown and still have 
wavelength two. The only way to describe it is to think of it as an upper cutoff 
mode that has a DC offset, which may also explain the transition of brighton to 
darkon. Why this should be so is not understood even slightly at this time. 

It remains to be seen whether an analogous structure exists in the softening 
lower cutoff mode; in fact, there may exist many more solitons in the "simple" cutoff 
modes. This complexity in the simplest of modes, which can be described by a 
relatively simple evolution equation in which all but the second spatial derivatives 
can be ignored should lead one to expect a far more difficult time dealing with 
intermediate modes, the subject of the next chapter. 
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State of Lattice After Brighton Decay 
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Figure 111.24. Final darkon like state after brighton decay. 
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Lattice Position 



IV. THE NONLINEAR LATTICE II: INTERMEDIATE MODES 



A. THE "LAMBDA FOUR" MODES. 

The numerical study of solitons in intermediate modes requires first that the 
modes themselves be understood. It turns out, at least for the "lambda four" modes 
(those whose wavelength is four), this is not a trivial matter. The plural is used 
because there exist a variety of manifestations, each of which has its own 
characteristics and preferred region in the drive plane. The lambda four case was 
chosen to follow the cutoff modes in order of study because it had been much 
observed experimentally. The mode has a natural linear frequency given by 



There are three known stable states for the lambda four pure mode, shown in 
Figures IV. 1, IV.2, and IV.3. During the study of these modes, they each acquired 
a descriptive vernacular name. The first to be studied was the "plus zero minus 
zero" or " + 0-0" mode, shown in Figure IV. 1. It was during the study of this mode 
numerically, and in parallel on the experimental lattice, that the preferential 
existence of the "plus plus minus minus" or "+ + -" mode in certain regions of the 
drive plane was discovered. The third stable version (there are presumably many 
more...) was discovered during the numerical study of the " + 0-0" mode drive plane. 




IV. A. 1 
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Fig IV. 1. A lambda four mode -- the "+0-0" Mode. 
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Fig IV .2. A lambda four mode -- the Mode. 
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The drive plane for the " + 0-0" state of the lambda four mode was studied 
extensively as a prelude to the study of kink drive planes (the drive planes for the 
other two states have not been studied systematically yet). Figure IV.4 shows the 
drive plane study results for one set of initial conditions. The lattice used was of 
forty elements. It is probable that, for a lattice of different size, the exact positions 
of the lines shown may vary, but how or whether they actually vary is unknown. 
The very linear nature of several of the boundaries of the drive plane region of 
stability was quite a surprise; the reason for the linearity is not understood. 

The behavior of the lattice at each of the boundaries of the drive plane region 
of stability proved to be very interesting, and provided much evidence to aid in the 
understanding of the lattice’s characteristics in general that proved valuable when we 
moved on to the study of kinks. 

The behavior of the lattice at the lower boundaries proved to be the most 
complex. For all drive frequencies less than 0.98, the lattice went unstable when it 
reached the boundary shown. The instability was slow to develop. This instability 
is caused by the commencement, when amplitude is lowered to the curve, of very 
slow and growing motion of the line of nodes (the 0’s in " + 0-0"). This is also the 
change that characterizes lattice behavior when the lower limits of the region of 
stability of the uniform mode are reached for drive frequencies greater than 0.98; 
however, the instability does not grow without bound there and will be discussed 
separately. 
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Drive Plane for the +0-0 Mode 




Fig. IV.4. Drive plane for " + 0-0" Mode. 
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Angular Frequency of Drive 



That this in phase motion of nodes begins to occur and that it grows can be 
understood qualitatively by considering the potential energy in which the elements 
rest. In the region of stability, the nodes are stable because they rest at the bottom 
of a potential energy well caused by the coupling to the adjacent high amplitude 
elements. In the pure mode, these adjacent elements have identical amplitudes, so 
that the two opposite coupling forces exactly cancel. If we consider a slight 
perturbation of the position of one of these nodes, it will feel less force from the 
element in whose direction it moved, since the relative displacement has decreased. 
Conversely, it feels a greater force from the other element. The net force acting on 
the node element when it is displaced differentially from its rest position is then 
restoring, and the nodes remain stable. In fact, the "numerical temperature" of the 
lattice could be determined by zooming in the amplitude scale until the nodes were 
seen to move -- a zoom of about 2 40 (or 10 12 , the double precision limit!) The 
motion is random, and corresponds by analogy to the thermal motion of crystal 
elements in their lattice positions. Now, as the drive amplitude is lowered, the 
parametric excitation of the nodes decreases. At the same time, however, the 
potential energy well in which they rest becomes shorter, since the amplitude of the 
antinodal elements is decreasing. Proceeding by analogy, considering the natural 
frequency of the node if it were to oscillate in the potential energy well as a simple 
harmonic oscillator, this change in well height would lead to a decrease in natural 
frequency. Using equation (I.C.13) for the threshold condition of parametric drive, 
one sees that a decrease in natural frequency leads to a quadratic lowering of the 
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drive threshold. While not rigorous, this discussion suggests how the lowering of the 
well leads to a lowering of the drive threshold at a rate faster than that at which the 
drive is lowered. Thus, at some drive amplitude, the threshold condition is met and 
node growth occurs. 

The reason the node growth continues until it destroys the mode (since this is 
a softening system, this merely means the lattice has at least one site which goes 
over the potential energy hump) is not well understood, although it may be related 
to the Benjamin-Feir instability discussed in Chapter II for cutoff lattices. When the 
nodes grow sufficiently to become commensurate with the amplitude of the antinodal 
elements, they continue to act initially as a unified group. Thinking of them as a sort 
of lower cutoff lattice embedded in the upper cutoff antinodal lattice, it can be seen 
that, for amplitudes greater than some threshold which depends on coupling and 
lattice size, the "nodal lattice" will become unstable to sidebands. For higher drive 
frequencies, the response amplitude is lower, so there should be a threshold 
frequency above which the Benjamin-Feir instability threshold is not reached. This 
appears to be what happens at drive frequency of 0.98. 

Above this threshold frequency, the nodes grow until a steady state is reached 
such as that shown in Figure IV.5. In a phenomenon similar to the weak instability 
which occurred for the N=2 lattice below the Benjamin-Feir threshold, a weakly 
unstable path to chaos exists at the second threshold line for higher frequencies (see 
Figure IV.4). This threshold occurs for decreasing drive amplitude , which is sensible, 
since it is associated with increasing parametric excitation due to a still decreasing 
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potential energy well from the antinotlal elements. Recall that, for the two element 
lattice, a threshold of weak instability was crossed which corresponded to the 
exceeding of a parametric drive threshold; in this case, the parametric drive 
threshold is for decreasing amplitude, and corresponds to the point where node 
motion begins. The chaotic state of the two element lattice occurred as one moved 
further from the drive threshold, as one crossed an apparently fractal boundary 
between the switching and non-switching basins of attraction; here, the switching 
and non-switching analogy corresponds to the switching between nodal lattice being 
driven by and driving the antinodal lattice. Although this connection between 
switching and non-switching is conjectural to some degree, since it is not absolutely 
clear exactly what is happening at the chaotic threshold, it is clear that as we move 
further from the threshold of parametric instability, a chaotic band begins at a 
definite threshold. 

Throughout this small chaotic region, the nodal lattice continues to maintain 
its identity and oscillates at a mean amplitude well below the amplitude of the 
antinodal lattice. Its motion is jerky and aperiodic, however, which typifies chaotic 
behavior. As in the case of the two element lattice’s approach to chaos, the 
development of the nodes as they approach the chaotic threshold can be illustrated 
well using frequency and phase domain visualizations. In Figure IV.6, the spectra 
of a nodal element are shown for drive amplitudes 0.075 and 0.056. The first is well 
above the threshold for node motion (which is 0.056); the spectrum is of thermal 
noise in the rough vicinity of the drive frequency. The second is taken early in the 
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growth of the nodes. The frequencies present are 1.10 and 0.87, which do not match 
the drive frequency of 1.0434, but are the frequencies about which the thermal noise 
centered in the first spectrum. The fact that the main frequency is higher than the 
drive frequency is consistent with the potential energy well idea, since this additional 
potential energy will increase the frequency of the nodal elements. The growth of 
the steady state oscillations of the nodal lattice as one lowers drive amplitude toward 
the chaos threshold can be seen in Figure IV.7, which shows phase portraits of 
several nodal and antinodal elements in their steady state at drive amplitudes of 
0.053 and 0.051. The trend, which continues as one lowers drive amplitude, is for 
the inner ellipses, which are the nodal elements, to grow outward, with the angle of 
inclination of the ellipses unchanged. At the same time, the outer ellipses, which are 
the antinodal elements (note they appear in opposite quadrants, since the antinodes 
are in an "upper cutoff" arrangement), fatten up, with their inner edges broadening 
and moving inward. As in the two element lattice case, this change is due to the 
greater amount of energy transferred from antinode to node during each cycle; there 
is as yet no switching of roles. 

When the drive amplitude reaches 0.045, the ellipses are just about to touch, 
as shown in Figure IV.8; the ellipses stretched a little further from this state and 
then chaotic behavior began, as seen in the lower half of Figure IV.8. In a manner 
similar to that of the two element lattice, the ellipses show large scatter. The spectra 
of this chaotic state are shown in Figure IV.9 for a node and an antinode. The 
characteristic spread spectrum of chaotic motion is readily apparent. However, as 
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Fig IV.6. Spectra for nodal element, (a) Thermal noise, (b) Onset of growth. 
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Figures IV.9 and IV. 10 show, this chaos turned out to be transient in nature 
(although it remained for many periods of oscillation). Figures IV. 10 and IV. 11 are 
the same as Figures IV.8 and IV.9, only they were taken after the transient chaos 
died out. Figure IV. 12 shows the continuation of the effects seen before when 
lowering drive amplitude; there were no transient chaotic states. Finally, at drive 
amplitude of 0.042, steady state chaos set in. This is seen in Figure IV. 13. Now 
there is clear penetration of the outer ellipses by the point of the inner ellipses, 
which I believe represents graphically the point at which switching of the nodes from 
driven to driving state occurs. 

An even more vivid graphical example of the route to chaos is seen in Figure 
IV. 14, which is a similar succession of phase portraits at various drive amplitudes, 
with the drive frequency set at 1.054 (higher than in the previous example, which is 
evident also from the lower threshold amplitude for chaotic motion). In the third 
portrait of Figure IV. 14, the lower point of the inner ellipse has clearly pierced the 
region of the left outer ellipse. Figure IV. 15 shows the same elements a few cycles 
later, with the idea of chaos very obviously portrayed. One point worth making 
about the change in shape of the outer ellipses is that, for all frequencies, when 
piercing occurs, the tails of the outer ellipses stretch out and bend down toward the 
coordinate (x) axis; this is because, as switching now occurs, these elements drive the 
nodal elements until they are nearly at rest in the equilibrium position (i.e., at the 
origin). This state was steady state (that is, it remained chaotic, and did not settle 
out into an ordered system); no transient chaotic states were observed at a drive 
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Fig. IV.9. Spectra for node and antinode showing transient chaos. 
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frequency of 1.054. 

When drive amplitude was lowered still further, the lattice made a transition 
from chaotic states to linear combinations of normal modes at about amplitude 
0.035. This corresponds to the regime wherein the effects of the nonlinearity in the 
equation of motion becomes insignificant. At somewhat lower amplitudes still, about 
0.025, the lattice motion died out, as the parametric threshold for steady state 
excitation was no longer exceeded. 

The upper right line in Figure IV.4 forms the boundary of the region of 
stability defined by a similar instability to that encountered at low drive amplitudes. 
When this line is reached, the nodes begin to move, but they are exactly 180 degrees 
out of phase with each of their node neighbors. The mechanism for this onset of 
node motion, while probably related to that which drives the in phase motion just 
described, is not understood. In all cases, the motion grows continually until the 
nodes and antinodes are commensurate in amplitude. Depending on position on the 
line of instability, various types of final states occur. For drive frequencies greater 
than about 1.04, the state that tends to arise is the ++- mode (Figure IV.2), 
although, with just a slight push beyond the line of instability, one can achieve the 
mode shown in Figure IV.3, which may be preferred for certain drive parameters 
(this has not been checked). For drive frequencies less than 1.04, the transient 
caused when the nodes become commensurate with the antinodes invariably causes 
the lattice to blow up due to the usual potential energy problem. What is most 
remarkable about the out of phase node motion boundary is that it is very linear all 
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Fig IV.ll. Spectra for node and antinode, after transient chaos. 
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Fig IV. 13. Steady state chaos -- spectrum and drive portrait. 
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Fig IV. 15. Phase portrait of steady state chaos at omega = 1.054. 
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the way to drive frequencies equal to the upper cutoff frequency, where it bends 
upward slightly. Why it should be so is one of the many aspects of the lambda four 
modes that is not understood currently. 

The final boundary of the stable + 0-0 region in the drive plane, the upper left 
curve, is characterized by a totally unexpected phenomenon -- isochronous motion, 
with respect to the drive, and DC excitation. When this boundary is reached, the 
nodes move rapidly to a fixed DC offset, which is matched by a DC offset in the 
antinodes, and stays there. If amplitude is further increased, the offset magnitude 
increases, but at the second curve, the antinodes become unstable, breaking up into 
complex modulated states which quickly lead to the potential energy problem and 
lattice breakup. The nodes and antinodes, when they move to the isochronous mode, 
may do so with an imposed modulation, depending on initial conditions as the 
boundary is approached, the modulation stays fixed once the steady state offset is 
achieved. Examples of modulated and unmodulated isochronously excited lattices 
are given in Figures IV.16 and IV.17. 

As is the upper right boundary, this boundary is very straight for a large range of 
values. In this case, for frequencies less than about 0.91, the curve bends upward; 
also as was the case for the upper right boundary, the reason for this profile is not 
understood. 

As this detailed look at the drive plane region of stability of the plain +0-0 
mode makes plain, the lambda four modes themselves, without solitons or other 
complications, exhibit quite complex and interesting behavior that is not fully 
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Fig IV. 16. Unmodulated isochronously excited lattice. 
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Fig IV. 17. Modulated isochronously excited lattice. 
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understood. The drive planes of the other forms of the lambda four lattice have not 
been studied, either and may present an additional range of interesting phenomena. 



105 



B. KINKS IN THE LAMBDA FOUR MODES. 



As one might expect, considering the complexity and number of lambda 
four pure modes, there exists a wide variety of kinks in these modes. Considering 
just the +0-0 mode, we can argue on purely geometrical grounds for the existence 
of four types of kinks, corresponding in effect to phase shifts of 90,180, 270 and 360 
degrees. The profiles they must exhibit, shown in Figure IV. 18, are based on the 
requirement that the entire lattice oscillate at the same frequency. The classification 
numbers given in Figure IV. 18 will be used throughout the following to clearly 
identify the various kinks; the profiles are based on a softening lattice, which was the 
case considered in the pure mode drive plane analysis described in the previous 
section. The main emphasis in this section is indeed on softening +0-0 kinks, 
although results are also presented for hardening kinks and kinks in the other 
lambda four modes. 

All of the kink types shown in Figure IV. 18 have been shown numerically to 
exist. Representative examples are shown in Figures IV. 19 through IV .22. 

In addition, Figure IV.23 shows a Type LA kink, which is an antisymmetric version 
of the Type I shown in Figure IV. 18. Figure IV.24 demonstrates clearly the 
distinction between symmetric and antisymmetric kinks. Presumably, antisymmetric 
versions exist of the other kinktypes as well, although this has not been verified. 
They are seen to be in good agreement with the expected general profiles. 
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7 ig. IV.19. Numerical results for a Type I Softening +0-0 Kink. 
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Type II Kink (+0-0 Mode) 
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Fig. IV.20. Numerical results for a Type II Softening +0-0 Kink. 
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Fig. IV.21. Numerical results for a Type III Softening +0-0 Kink. 
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Type IV Kink of +0-0 Mode 




Fig. IV.22. Numerical results for Type IV Softening +0-0 Kink. 
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Type IA Softening +0-0 Kink 
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Fig. IV.23. Numerical results for Type IA Softening +0-0 Kink. 
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Fig. IVJ24. Symmetric and antisymmetric +0-0 kinks. 
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In order to make a clear comparison between kink behavior and normal mode 
behavior, a drive plane study was conducted of a Type I kink. Figure IV.25 shows 
the results of this study, with the outline of the +0-0 pure mode drive plane overlaid 
for comparison. Several striking results were obtained in the course of this drive 
plane exploration, as expected, given the complexity of the modal drive plane 
boundaries. In particular, the presence of the kink produced dramatic, localized 
effects which appeared in many cases to compete with the effects in the wings, which 
generally were identical to the ones discussed in the pure mode drive plane results 
of the previous section. This is sensible, since in the wings, that is, far from the kink, 
the lattice should be expected to behave as a pure mode. 

Looking first at the lower boundary of the stable region for the Type I kink, 
we see that the general behavior is similar for drive frequencies above 0.87. In 
particular, however, we note that the frequency at which node growth stabilizes has 
increased from 0.98 to 1.02. This is apparently due to the action of the kink in 
driving the node growth, which occurs most rapidly in the vicinity of the kink. 
Chaotic motion begins at roughly the same spot, with drive frequency of 1.06 and 
amplitude of 0.041. A particularly striking phenomenon that was observed was the 
migration of the kink to the right side of the lattice (from the middle) when the 
drive amplitude was lowered to 0.037, and the chaotic motion was strong. Figure 
IV .26 shows this kink after its migration. 

Far to the left, below drive frequencies of 0.87, however, the boundary bends 
upwards form that observed for the pure mode. The reason for this was not initially 
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Fig IV.25. Drive plane for Softening Type I +0-0 Kink and Pure Mode. 
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Fig IV.26. Test kink after migration. 



116 



clear, until suddenly at drive frequency 0.795, the kink went through a severe 
transient in which it made a transition from symmetric to antisymmetric states, with 
the latter being much sharper, that is, occupying two vice four lattice sites. It 
became clear then that the reason for the raised boundary is that, as we lower drive 
frequency and increase response amplitude (since it is a softening system), the kink 
tends to sharpen up. By analogy to the NLS kinks for which theoretical descriptions 
are well established, this sharpening was expected; what was unexpected was the 
sharp transient between kink geometries. This being the case, it appears that the 
symmetric kink can survive even at very low drive frequencies, but only at elevated 
drive amplitudes (which of course precludes node motion). If the amplitude is 
lowered, the nodes begin to move so as to bring the kink to its more stable (or lower 
energy) state, the antisymmetric state. As a qualitative check for this, it was verified 
that the antisymmetric kink was stable much closer to the original boundary of pure 
mode instability. 

For drive frequencies above 1.06, which is close to the linear frequency of the 
upper cutoff mode, the line corresponding to out of phase node growth follows 
closely that of the pure mode, including the bending upward in the vicinity of drive 
frequency 1.09, which is the linear upper cutoff frequency. However, below drive 
frequency 1.06, the kink’s stability boundary deviates sharply from the mode’s. From 
that point all the way left to the isochronous curve (which matched that of the linear 
mode very closely), the boundary amplitude increased in a stepwise fashion. This 
very unusual and inexplicable result was checked carefully at the points indicated. 
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although it is possible that the boundary exhibits curvature in between the data 
points shown. This deviation from the pure mode case extends to the type of 
behavior as well as its location. In the small region above drive amplitude 0.11 and 
between frequencies of 0.99 and 1.04, there exist stable states with out of phase node 
motion; a steady state example is given in Figure IV.27. Observation of the 
development of such a steady state makes clear what the mechanism of the 
departure from the pure mode is. The kink apparently drives the nodes closest to 
it, so that, for sufficient drive amplitude, the node motion is propagated down the 
line of nodes and eventually grows. The node motion is thus a radiative event, and 
the node motion is not precisely out of phase as it is when the "out of phase 
instability" boundary is reached (that is, the upper right boundary of the pure mode 
drive plane). When the drive amplitude is increased at drive frequency 1.04 to just 
under the original out of phase instability curve, the lattice does indeed transition to 
the + +-- mode in a complex (three kinks) state (Figure IV.28), suggesting that 
boundary still exists, but is just not reached before the kink drives the lattice out of 
its original steady state and into another. For drive frequencies below 0.99, the 
effect of the kink, which is now sharper and of greater amplitude, in driving the 
lattice is too great for a steady state to be reached with the nodes still small 
compared to the antinodes; the system "blows up". 

Not all of the kinks observed numerically were of the "pure" types shown in 
Figures IV.18 through IV.23. Frequently, after some perturbation or change in initial 
conditions the lattice would, after a transient of varying length, emerge with one or 
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Fig IV .27. Steady state upper cutoff node motion in Type I kink drive plane. 
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Fig IV.28. Three kinks in + +-- mode 
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more complex kinks. An example is shown in Figure IV.29. Using mismatch 
between left and right domains as our classification guide, we see that this is a Type 
II kink (that is, we find that antinodes in the right hand domain are where they 
should be if there had been no kink, but they are "up" when they should be "down" - 
- Type II kink). Denardo has suggested that these complex versions of simple kink 
types be called "excited kinks", which is an attractive nomenclature, since it recalls 
the particle like qualities of solitons and confers on these states the characteristics 
of a particle which is in an excited energy state. Another way of viewing these 
phenomena, which ties in with the discussion later in this chapter of domain walls, 
is as kinks with inclusions, which might be termed the solid state analogy -- viewing 
these solitons as analogous to transitions in crystal lattices. Both views suffer from 
a lack of mathematical underpinnings, but they offer complementary insights into 
what is happening in the lattice. 

In Chapter II, the symmetry properties of NLS solitons were noted and 
verified. It is a matter of great interest to determine the point on the dispersion 
curve, if there is one, about which this symmetry is based. Unfortunately, this 
question has not been resolved. However, some tantalizing hints have been found. 
In particular, it has been found that, if one takes a Type I softening kink of the +0-0 
mode and performs the following manipulations on it, it transforms into a positive 
energy hardening kink: 
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IV.B.2 
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where we have relied on the fact that we set 

« 0 -l. IV.B.3 

The hardening kink thus created has a span that is similar to that of the original 
softening kink, as seen in Figure IV.30. In the cutoff modes, such a transformation 
would convert from breather to kink, but here it only takes us from a negative 
energy to a positive energy kink - the wings here still have finite amplitude. Thus 
the relation of this phenomenon to the symmetry between breathers and kinks is not 
confirmed; in fact, no attempt has been made yet to determine whether lambda four 
breathers even exist. 

A classification of hardening +0-0 kinks can undoubtedly be made along the 
lines indicated above for the softening case; here we have focused almost exclusively 
on the softening case in order to take advantage of the experimental work going on 
with the softening pendulum lattice. An additional type of hardening + 0-0 kink has 
been observed, and is shown in Figure IV.31. 

One phenomenon that was observed in a hardening lattice that has not been 
observed elsewhere yet, and which has tremendous possibilities, was a temporally 
phase shifted kink (Figure IV.32). Whereas all of the previous nonlinear structures 
studied in this work and its predecessors have relied on spatial phase shifts and 
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Complex Type II Kink in +0-0 Mode 
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Fig IV.29. Complex Type II +0-0 kink. 



123 



Lattice Position 




(qjB) eprundwv 



c 

o 



U 5 

o 

CL 

0 

O 

's 



Fig IV .30. Example of hardening +0-0 kink. 
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Fig IV .32. Additional type of hardening +0-0 Kink. 
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Hardening +0-0 Kink with Phase Shift 

(T wo middle kink elements lag by 8 deg) 
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Fig IV.33. Hardening +0*0 kink which relies on Temporal phase shift. 



126 



Lattice Position 



amplitude modulation to maintain monofrequency response, which is the bottom line 
requirement for steady state to exist, this kink relied additionally on an eight degree 
temporal phase lead in one of the elements to maintain monofrequency response. 
There is no reason to assume that there isn’t an entire range of phenomena relying 
on similar phase shifts, but thus far only the one example has been observed (to my 
knowledge). 

The extreme richness of behavior encountered in the +0-0 mode led to a 
conscious decision to focus on it to the virtual exclusion of the other lambda four 
modes. However, there is every reason to expect a similar degree of complexity to 
govern the behavior of these modes as well. The large number of additional modes 
intermediate between the cutoff modes have been studied almost not at all, although 
many kinks have been seen in them as a by product of the work presented here. 
Specifically, kinks have been observed in modes with wavelengths of three, five, six, 
seven, and eleven. The phenomenon of kinks appears to be ubiquitous. 

C. DOMAIN WALLS IN THE NONLINEAR LATTICE. 

The phenomenon of domain walls is well known from the study of crystal 
structure, electromagnetism, and superconductivity. A domain wall in a lattice is a 
(uaually sharp) boundary between two different domains. An example would be a 
boundary between the upper cutoff mode and the lambda four mode. Typically, the 
two domains are independent of each other and of the wall beyond a very short 
interaction distance; that is, domain walls are local phenomena. They have been 
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observed experimentally in the pendulum lattice, although systematic study only 
began recently, in parallel with this work. There is very little theoretical 
understanding of the phenomenon, however, especially when the transition between 
domains takes place over a finite length scale. 

While domain walls had been a topic intended for research at the end of this 
work, if at all, events intervened when they appeared on their own during an 
investigation of lambda four kinks. When the positive energy (Type II) lambda four 
kink was first demonstrated, the form obtained was not that given in the previous 
section. In fact, it was an "excited" state, or a state with inclusions, depending on 
how one chooses to view the phenomenon. Figure IV .33 shows the original state, 
and Figure IV.34 shows what appeared when the indicated elements were removed 
(by dumping the data to a file and then editing the file). As guessed, the removal 
of the indicated elements did not adversely affect the stability or identity of the kink. 
The remarkable flatness of the "kink", as it was initially viewed, suggested strongly 
that it might instead be a domain wall. To test this hypothesis, the elements 
constituting the flat region of the structure were removed and the lattice restarted. 
Not too surprisingly, a stable positive energy kink of normal size resulted. However, 
it was also possible to extend the flat region by placing identical elements (i.e., 
elements with same amplitude and velocity) in the middle of it, and then removing 
the original primary domain. This resulted in a stable ne gative energy kink of the 
upper cutoff mode! 
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Fig IV.33. Lambda Four Type II Kink with inclusions. 
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Fig IV.34. Same kink as Figure 34, with indicated elements removed. 
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The simple elegance of this result strongly suggests that domain walls might 
profitably be viewed as regions where a negative energy kink in one domain and a 
positive energy kink in the other exactly match some set of boundary conditions, 
such that the pair of kinks is stable, and the domains on each side are unaffected by 
the presence of the other. It is interesting to realize that the use of periodic 
boundary conditions in my model made this conclusion more evident, since one 
always needs and even number of domain walls in order to meet the periodic 
boundary conditions. Thus the idea of either domain as an inclusion in a kink of the 
other is visually evident, whereas if domain walls had been studied as structures in 
their own right from the beginning, the connection might not have become apparent. 

Since that auspicious beginning, the only significant result obtained concerning 
domain walls (excepting the reaction to media nonuniformities, which will be 
discussed in Chapter V), is that they too are ubiquitous. This fact should not be 
surprising, since kinks, the building blocks (evidently) of domain walls, are 
ubiquitous. In fact, the wide variety of domain walls observed to date, including the 
extreme case of upper cutoff/lower cutoff domain walls, suggests that, for any 
positive energy kink, there will exist a domain wall solution with any mode that 
exhibits a negative energy kink, whose amplitude would be higher than the original 
mode, and vice versa. So, for example, in the hardening case, one would expect an 
upper cutoff positive energy kink would have domain wall solutions with any other 
mode which the lattice can support (i.e., limited only by lattice size and mode 
stability constraints). 
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The further study of intermediate modes is obviously necessary, and domain 
walls are one of the most important areas. Their potential bearing on many 
technologically important areas in solid state physics and critical point phenomena 
demands a close examination. Moreover, the need for a comprehensive theoretical 
treatment that allows all of the phenomena observed to date is critical. 
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V. EXTENSION OF THE BASIC RESULTS. 



A. THE TWO DIMENSIONAL LATTICE. 

As a natural extension of the work reported to this point, a two dimensional 
model using the same physical ideas was developed. Only some very preliminary 
work has been done, but enough has been discovered to warrant future work in this 
area. While no theoretical work has been done for the two dimensional case yet, it 
will be useful to speculate a little on the likely course the theory will follow, since 
it allows an interpretation to be made of the results presented that at least seems 
reasonable. 

The model used still featured nearest neighbor interactions only, but now in 
two dimensions. Diagonal interactions are specifically ignored, however. The exact 
equation of motion is given by 

V.A 

Here m and n represent the row and column number, respectively. We can see 
quickly, by analogy to the one dimensional case, that the "upper cutoff" case, where 
each element is exactly 180 ° out of phase with each of its four neighbors, has a 
linear frequency given by 
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It is also clear that, in the case where the lattice is upper cutoff along one of the 
axes of the lattice and lower cutoff along the other, that 
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V.A.3 



which of course is the linear upper cutoff frequency for the one dimensional lattice. 
This is obvious, since the lattice sees no coupling in the lower cutoff axis’ direction. 
In fact, the elimination of diagonal interactions effectively decouples the two 
orthogonal directions x and y, so that the linear dispersion relation can be written 
by direct analogy to that of the one dimensional lattice: 
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Considering the cutoff cases, of which there are now four, we see that, again 
by analogy to previous work, 
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where m and n are either 0 or 1. Without for the time being proceeding to formal 
proofs, we note the similarity of this equation with the one dimensional case, and 
speculate that the system may be represented by a two dimensional NLS equation. 
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Accordingly, it is reasonable to suppose that one might find soli tons in the two 
dimensional cutoff lattice resembling the NLS solitons observed in one dimension. 
It may be that the soliton is truly two dimensional, or it may be coupled to a pure 
mode in the orthogonal direction. Figure V.l shows some of the structures one 
might hope to find in the cutoff two dimensional lattice. 

Not surprisingly, the two dimensional lattice displayed fascinating behavior 
from the very beginning. The initial conditions used for most of the earliest work 
consisted of sinusoidal modulations in the x and y directions, with exactly one 
wavelength spanning the x and y directions. In all of these cases, the lattice 
displayed behavior that could be divided neatly into three time scales. On the first 
time scale, of a hundred or so periods, consisted of a rapid disordering of the lattice 
as the initial disturbance radiated energy (now in two dimensions, so the radiation 
is more complex). Over the next several hundred periods, in every case the lattice 
resolved itself into a small number of domains of the original modulated cutoff 
mode, with each of the domains being 180° out of phase with its neighbors. Figure 
V.2 shows the end of such a period for one set of initial conditions. On the final 
time scale, a mechanism which is not well understodd but acts in a fashion similar 
to surface tension (see below) did indeed reduce the lattice to a single cutoff 
domain. This time scale was often very long, since in some cases the surface tension, 
or net difference in force sensed by each domain, was often quite small. Figure V.3 
shows the same lattice as seen in Figure V.2, late in this evolution. 
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Fig V.2. A 2D lattice at the end of the second time scale. 
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Fig. V.3. The same lattice, late in the final time scale. 



138 



The results of these initial investigations made it clear that, if kinks were to 
exist, they would have to exist in some geometry which gave neither domain (there 
are effectively only two, since they are 180° out of phase) a dominant position. A 
few obvious geometries are shown in Figure V.4; all of them did indeed support 
kinks. Figure V.5 shows a very striking example, which resembles two NLS negative 
energy kinks, one in the x and one in the y directions. This result was very 
encouraging, if not unexpected, for it strongly supports the notion that the theoretical 
description of the two dimensional lattice will be very similar to that of the one 
dimensional lattice. Of course, this rosy picture would immediately fall apart if 
diagonal interactions were allowed in the model, as cross derivatives would then 
abound in the equations of motion! 

Figure V.6 shows the same kink as seen in Figure V.5, but with a complex 
structure that resides at each of the intersections between the x and y kinks. This 
was actually the first two dimensional kink seen, and was used as the starting point 
in getting to Figure V.5. The structures seen are very stable, and they are certainly 
not yet understood. 

The final quick foray into the two dimensional arena which was undertaken 
was an effort to determine if there existed stable breathers. Initially, it was felt that 
such structures might not exist, due to the existence of diffraction effects. However, 
they were found to exist; in fact, it turns out that if one starts with a lattice 
completely at rest and then kicks one element (preferably the middle one for 
visualization purposes) with a reasonable amplitude (i.e., enough that nonlinearities 
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Fig. V.4. Possible geometries for stable 2D kinks, where surface tension is zero. 
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Fig. V.5. An actual 2D lattice with symmetric kinks. 
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7 ig. V.6. The same lattice with complex structures at the kink intersections. 
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can play a significant role), and in fact a great deal of radiation occurs initially which 
leaves behind a very clean two dimensional breather, seen in Figure V.7. It is 
actually easier to get a two dimensional breather than a one dimensional one, 
because the radiation spreads out in two dimensions and is much more rapidly 
attenuated. Even in the free case, a very clean breather can be obtained, since the 
energy radiated away as the breather seeks its equilibrium shape has to spread itself 
over the entire two dimensional span of the lattice, resulting in very low amplitude 
radiation passing through the breather at any given time (although of course there 
is always some). 

A theoretical consideration which was suggested for the first time by the results 
obtained with the two dimensional lattice is the idea of "surface tension" at domain 
boundaries. Suppose for the moment that two similar domains are connected by a 
kink of some two dimensional shape. An example which will motivate the discussion 
is given in Figure V.8. The question arises whether one or the other of the domains 
will dominate the other; that is, will one domain gradually force all of the elements 
of the other to shift modes? As it turns out, this is indeed the case. In fact, for two 
domains of similar amplitude (i.e., identical domains separated by kinks), the domain 
which is concave will dominate the convex domain, regardless of the size of the 
domain. 

That this is the case is reasonable, when the issue is considered as analogous 
to the phenomenon of surface tension. Since the concave domain exerts more force 
per boundary element on the convex domain than vice versa (because it has more 
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Fig. V.8 Surface tension in two dimensional domain boundaries. 
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elements interacting at the boundary), there is an effective surface tension, or net 
force, acting on the convex domain. Since, unlike the case of an air filled bubble 
being compressed in water, there is no reaction force which builds up as compression 
takes place, the convex domain is gradually "compressed" away. It should be noted 
that this phenomenon depends on the fact that coupling of an element to its home 
domain is negligible, so that one can simply add up the coupling forces to the 
opposite domain to determine the dominant domain. 

Presumably a modified version of this rule would be operative when dissimilar 
domains are in contact; however, it is not immediately clear what exact form the rule 
would take after the differences in mean amplitudes are accounted for. In any case, 
no numerical data for such cases has so far been collected, so it remains a truly open 
question (as does virtually everything having to do with two dimensional lattices 
given by (V.A.l)!). 

This brief extension of the lattice model into two dimensions leaves much left 
to explore. It is clear, though, that there are many new phenomena likely to reward 
such efforts. The model also permits study of surfaces, including radiation and finite 
structures, providing that sufficient computer power is available to provide a small 
enough discretization of the surface (that is, more memory is needed, since the 
model is limited on the machine used to a 40x40 lattice, which is not adequate to use 
for a surface model). 
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B. THE NONUNIFORM LATTICE. 



As a second attempt to extend the work reported in Chapters III and IV, an 
investigation was made into the effects of nonuniformities on lattice behavior. To 
some extent, this work was motivated by results obtained on the experimental lattice, 
since it was assumed that that lattice was nonuniform, and the effects and indeed 
types of nonuniformities needed to be understoodif the experimental work is to be 
joined to the numerical. Additionally, it was desired to test a method of removing 
the effects of nonuniformities which had been tried by the group at UCLA with 
which we closely worked. This method consisted of first measuring the amplitudes 
of each element in the lattice in a pure cutoff mode, then dividing these results into 
the measured lattice amplitudes in the presence of solitons. While this method was 
based on intuition, it proved to give smooth results from results that had been 
difficult to interpret. This method of division greatly facilitated their study of 
solitons, but it needed to be verified before full confidence was placed in it. 

Originally, variations in coupling, natural frequency, and nonlinear coefficient 
were tried. The nonlinear coefficient was found to affect results only slightly. 
Accordingly, only results from nonuniformities in coupling and natural frequency are 
discussed in this chapter. Of these, natural frequency was the most studied, although 
only because of time constraints. 

The parameters were varied in three different methods - randomly, with a 
linear gradient, and with sinusoidal gradients. The program allowed the user to 
choose the type of variation and the amount. Also, it was during these investigations 
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that a waterfall display was added to the program so that a parameter could be 
instantaneously changed by one of the three methods just mentioned, and the impact 
viewed in a format which made changes easier to identify. 

Finally, before proceeding to discuss actual results, it should be noted that no 
attempt has yet been made to fit these and the experimental results into a theoretical 
framework. This work follows in the time honored tradition of forging ahead with 
experimentation in order to provide grist for the theoreticians’ mill, although I 
frankly would have preferred to have tended that mill myself, had time allowed! 

The first result was obtained early in this thesis research, using an older version 
of the program in which the nonlinear coefficient was allowed to vary. It was found 
that the division method developed by the UCLA group [Putterman, unpublished, 
1990] was valid. An example is shown in Figures V.9 through V.ll. The rough 
soliton-like structure in Figure V.9 was divided, element by element, by the 
amplitudes of the cutoff mode in Figure V.10 to produce the smooth soliton seen in 
Figure V.ll. The method was found to be valid for both coupling and natural 
frequency variations. One additional item noted, which served as the starting point 
for later work on nonuniformities, was that, for softening nonlinearities, the solitons 
invariably moved to the local minimum in the varied parameter. Since reduction in 
coupling leads to a reduction in frequency, this result means that, in general, 
softening kinks move to the local minimum frequency in the lattice. 

In order to characterize the behavior of kinks, breathers, and domain walls in 
various types of nonuniform lattices, we returned to this line of inquiry recently. 
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Fig. V.9. Soliton in a random coupling lattice. 
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Fig. V.10. Actual variation in coupling of lattice in Fig. 1, with pure mode shown. 
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Fig. V.ll. Same lattice, after the UCLA division method was used. 
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The initial results have been consistent, but as many questions have been raised as 
have been answered. The first result in the second round of nonuniformity work 
complemented the last result of the first round -- hardening kinks move to the local 
maximum frequency in the lattice. In the case of continuous gradients in either 
natural frequency or coupling, the hardening kinks moved steadily to the right 
(increasing frequency direction) until they reached the end of the lattice, where there 
was a discontinuous drop in frequency due to the periodic boundary conditions. 
Since there were always two hardening kinks, this provided a convenient way to 
collide kinks, since the leftmost kink inevitably caught up with the pinned right hand 
kink. As expected, these encounters produced annihilations of both kinks, leaving 
an undisturbed nonuniform lattice mode. 

In order to better study collisions, and to discern whether the motion was a 
diffusion process or a process of motion in a potential field, a sinusoidal variation 
of coupling was introduced, with one kink on either side of one of the peaks in the 
variation. In particular, we wanted to find out whether, when a kink arrived at a 
maximum frequency, it overshot the mark and exhibited oscillatory behavior. In 
fact, this did not happen at all, so we were also unable to observe two kinks passing 
through one another in this way (since they both stopped at the local maximum). 
Nevertheless, this sinusoidal parameter variation feature was found to be useful for 
separating kinks, since, if one places a minimum between two kinks, they will move 
in opposite directions. Since the variations could be turned off as easily as on, one 
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could separate two kinks using a sinusoidal variation and then turn off the variation 
to observe the cleanly separated kinks in a uniform lattice. 

We turned next to domain walls, anxious to see whether they moved 
preferentially in one direction or the other. In fact, no matter how hard we pushed 
the lattice, we could not get domain walls to move. The UCLA group independently 
reported similar results. Since one possible mechanism for this pinning of domain 
walls was that negative energy and positive energy kinks might tend to move in 
opposite directions (viewing domain walls in this instance as paired kinks), I 
immediately set after positive energy kinks, to see if they moved. Unfortunately, 
they also appear to be pinned. Why this should be so is not clear, and the fact that 
it is so can be interpreted as an argument against the notion that kinks and domain 
walls are intimately connected phenomena. On the other hand, it may simply be 
that the departure from normal mode amplitude by individual lattice elements is the 
key cause of soliton motion in a nonuniform lattice. If this were the case, and it 
doesn’t seem unreasonable to suppose that it might be, domain walls and positive 
energy kinks alike would not move, or would move so slowly that the experiments 
performed to date would not have detected the motion, since the changes are small 
in these cases. Another way of expressing the idea is to note that, in the negative 
energy kinks considered here (only a small subset of those possible), the envelope 
one would draw of the lattice often changes sign across the kink, and always makes 
a major transition in the kink region. In the positive energy kinks and domain walls, 
and possibly "darkons", the envelope only changes slightly, and retains the same sign 
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throughout. One possible way of checking this idea is to measure the speed of 
motion of negative energy kinks as a function of amplitude, since presumably if the 
hypothesis were valid, they would move more rapidly as amplitude increased. In any 
case, this area is one which might fully occupy a future thesis student. 

For certain types of abrupt changes in lattice parameters, such as a strong 
discontinuous drop at the lattice ends, solitons were actually created when the 
parameter change was introduced. For example, in Figure V.12, we see a pair of 
kinks just prior to and just after the introduction of a sinusoidal modulation of 
natural frequency. When the modulation is introduced, a kink/antikink pair arises 
at element 90 and leftward (element 90 is a frequency maximum). Shortly after this, 
the leftmost kink of the created pair undergoes an annihilation event with the 
original kink, leaving behind a kink of similar type located where the original kink 
would eventually have moved in any case -- at the frequency maximum. Whether 
this is in effect the mechanism for all kink motion remains to be seen. 

Another way one might attempt to explain the difference between negative 
energy kinks, which move, and positive energy kinks and domain walls, which don’t, 
is to assert that very sharp solitons are pinned, whereas broad solitons are free to 
move. The ambiguity arises from the unfortunate reality that only a limited number 
of cases have been tried so far, and for all of the domain walls and positive energy 
kink cases, the structures are very sharp. It may be that, when a broad domain wall 
is obtained (if they exist), it will move exactly as the negative energy kinks do. 
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Fig V.12. Transition by kink creation and annihilation. 
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Finally, a quick look at breathers revealed that they do not seem to move at 
all, or even to be affected by the variations, except in this peak amplitude. This may 
again be due to sharpness effects; certainly the envelope of a breather is a rapidly 
changing function. This result is striking, for breather motion has been observed 
experimentally in water trough experiments (Wu, et al. [1984]). 

The best that can be said after this very preliminary look at the behavior of 
nonlinear lattices with nonuniformities is that the behavior promises to be as rich 
and varied as the uniform lattice, with many new phenomena undoubtedly awaiting 
later researchers. It would be very interesting to study the effects of nonuniformities 
in two dimensional lattices, but that’s another entire area of research.... 

C. THE TODA LATTICE. 

In some ways, the most important contribution this thesis makes to the study 
of dynamical systems is the introduction of a highly interactive and readily adaptable 
modeling tool for oscillating systems. The modeling should be limited to oscillating 
systems, since the errors incurred in using a simple Euler’s method derivative may 
accumulate in other types of physical models (Appendix A). The program offers the 
investigator quick access to an accurate model of his system of interest which allows 
him to view system dynamics in the time domain, frequency domain, and phase 
domain. Additionally, a full suite of file handling routines is provided so the 
investigator can save states of his system, or he can save sequences of time sampled 
data or spectral data in files while observing system behavior. 
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In order to demonstrate the utility and flexibility of the program, a model of 
the well known Toda lattice (Toda [1971]) was developed. It required approximately 
fifteen minutes to completely alter the program, converting it from a model of one 
physical system to another! At the conclusion of that time, I was taking data exactly 
equivalent to the types of data used throughout this thesis. The Toda lattice is a 
lattice where nearest neighbors interact according to an exponential interaction 
potential 



Here x n is the usual variable representing departure from equilibrium position of the 
nth element. 

The free system given by (V.C.2) suffers from the same problem that all free 
systems do in numerical work - it is often difficult to distinguish solitons from 
background radiation, which of course remains for all time due to the lack of 
damping. Accordingly, we desire to stabilize the structures of interest by considering 
damped and driven systems. In the case of the Toda lattice, there are many possible 
ways to add damping and drive. Geist and Lauterborn [1988] used a sinusoidal 
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which yields an equation of motion (Kuusela and Hietarinta [1990]) 
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driving force acting only on the first element of the lattice, and viscous damping. 
Since this is quite distinct from the parametric drives used in the lattice models we 
have studied so far, this method was not used. Instead, three different drive schemes 
were developed. The first two consisted simply of parametric drive via each of the 
two parameters of the system (a and b); the last was more elaborate. In this scheme, 
the lattice elements were each placed in an external potential well analogous to that 
provided by gravity in our previous lattices; the parametric drive was viz the natural 
frequency of the element resulting solely from its interaction with the potential well 
(in other words, the same parametric drive was used that was used in the other 
models of this thesis). 

Mathematically, the first two methods are given by 

x n -(a+2r)cos2<j3t)(e V.C.3 



and 
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The third method is given by 

j K n -a(e'^ Xn <>(J: "‘ , " x ' ,) )-((»)Q+2Ticx)s2a)r)jc n . V.C.6 

As before, we arbitrarily set the natural frequency to 1 in the program, without any 
loss of generality. 

All three of these methods were coded up and run on the computer. The first 
two have proven to be very difficult to control, with the slightest changes in drive 
or dissipation producing large results. It was possible to get some stable states, 
however, with one of them shown in Figure V.13. There are two lattices shown in 
the figure, because the stable state is, more accurately, a quasiperiodic state where 
the two states shown alternate back and forth, with a transient during the alteration 
which shows beating between the left and right hand sides of the lattice. This cyclic 
behavior is stable, hence the lattice could be termed quasistable. However, it seems 
to me that the main utility of the "pure" Toda lattice is its mathematical tractability 
(it is exactly integrable and has been studied extensively in the mathematical 
literature); it is difficult to imagine a physical system which obeys the equations of 
motion. 

More physically interesting is the third type of Toda lattice model, which can 
be interpreted physically as particles with exponential interaction potential resting 
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Fig. V.13. Two states of quasiperiodic Toda lattice with gamma drive. 
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each in its own external potential well. A physical realization of this is the solid 
state (this interpretation is due to Larraza). The exponential interaction potential 
is still questionable, but it is at least a feasible model. The third type of model was 
also the most practical numerically, in the time I was able to allot to the Toda lattice 
study. It should be noted that it may merely be a poor choice of parameters that 
resulted in the difficulty in using the first two models, and a more detailed study will 
undoubtedly prove useful. The first result obtained with the Toda lattice in an 
external potential was a negative energy kink similar to the cutoff kinks in the NLS 
lattice, in a travelling upper cutoff mode . That is, the mode, with its kink in train, 
moved from right to left continuously, at a fixed rate. This kink, could not be 
stabilized as a standing wave for system parameters close to the original parameters 
which yielded the results, although there is no reason to conclude that standing kinks 
might not exist elsewhere in the parameter space. In fact, for similar reasons we 
cannot really conclude that there are not travelling wave solitons in the nonlinear 
lattice model that the rest of this thesis addresses; we just haven’t observed then yet. 
Figure V.14 shows the profile of this travelling wave domain wall at one instant in 
time. It maintains this form as it moves to the left. 

The results shown in Figures V.13 and V.14 are only the tip of the iceberg, and 
are shown more to demonstrate the flexibility of the program and the computational 
ease with which a new line of inquiry can be undertaken. For example, it was found 
that lattices which were driven in the coupling term were very difficult to control 
and exhibited at best very limited stability. In any case, one of the beauties of this 
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Fig. V.14. Travelling domain wall in Toda lattice with external potential. 
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work is that the limiting factor in one’s progress is one’s understanding of the physics 
of lattices; there is no reason to lose productivity to the computer, since the 
modeling technique is well developed and very easy to use and mosify. It is to be 
hoped that a follow on thesis student will conduct a thorough study of all three Toda 
lattice models (and perhaps others), in conjunction with an attempt to understand 
the theory of Toda lattice solitons. There are some results in this area in the 
literature, but very little is known compared to the simpler lattice models. 
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VI. CONCLUSIONS AND PROSPECTS FOR FUTURE WORK 



It should be quite evident at this point that the behavior of nonlinear lattices 
is extremely complex. We have shown that kinks exist in virtually all modes of 
lattice vibration, and domain walls appear to exist between all compatible modes (by 
compatibility, we mean all modes where the higher amplitude mode has a negative 
energy kink and the lower amplitude mode has a positive energy kink). Neither of 
these two assertions has been rigorously proven, but the evidence for them is 
encouraging. Additionally, breathers in the cutoff modes were found to be very 
robust, and the remarkable self -focusing phenomenon discovered experimentally by 
Denardo [1990] was found numerically as well. Breathers in the Toda lattice and 
in the two dimensional analog of the basic lattice we discussed were also found with 
relative ease, further suggesting that breathers are common phenomena. 

What needs much more work, however, is the theory underlying these many 
nonlinear structures. The NLS theory for the cutoff modes is satisfactory, but it 
covers only the endpoints of the dispersion curve. It is hoped that, underlying this 
broad range of interrelated phenomena, there may exist a common equation of 
evolution which will prove more useful than the current favorite (the Korteweg-de 
Vries Equation), which has limited physical application. If such a general treatment 
that includes the NLS as a special case and which treats all intermediate modes can 
be found, it would find broad application in the physical sciences, since the 
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underlying equations of motion of the systems considered here are frequently 
encountered. Examples include such solid state topics as ferromagnetism, 
superconductivity, and nonlinear optics, as well as topics in plasma physics, 
cosmology, and particle physics. 

There are many questions raised during the course of this work which can be 
profitably addressed by follow on researchers using the model developed here; 
resolution of these questions will undoubtedly ease the task of achieving the 
theoretical understanding desired. It is still an open question whether domain walls 
and kinks are distinct fundamental structures or whether domain walls are always 
matched kinks (the latter is my view, but it is by no means unanimously accepted). 
In addition, there are several classification problems concerning kinks in 
intermediate modes, which may be useful to address as a means to guide the work 
of the theoreticians. And it is not yet known whether there exist any breathers in 
the intermediate modes. Indeed, the distinction between breathers and kinks, which 
was made during the study of cutoff mode solitons when the distinction was clear, 
may not be a useful distinction in the intermediate modes. And, it is unclear why 
the NLS theory works so well far outside the limits of its validity. 

This work focused mainly on kinks and domain walls, after the initial 
verification of breathers’ existence according to the NLS theory. It is left to the next 
student (who has already commenced work) to explore breathers in greater detail; 
it is hoped that the results of this later work will dovetail with those presented here 
to clear up the basic questions. 
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We only scratched the surface in this work in the study of two dimensional 
lattices and of lattices with different equations of motion (i.e., the Toda lattice). The 
work done was done primarily to flex the muscles of the computer program to verify 
the ease of adaptation which was one of the design goals of the model. The results 
that were obtained, however, were very interesting, and make the pursuit of these 
lines of inquiry highly desirable. In particular, the two dimensional case provides an 
opportunity to explore the effects of diffraction on soli tons. Just from our 
preliminary look, this effect appears to be very important, as it accounts for the 
unexpected stability of two dimensional breathers (which were thought before to be 
impossible to obtain). 

In conclusion, we were fortunate in this research effort to be able to make 
many interesting discoveries and to bring together conclusively some theoretical and 
experimental work. The results obtained have great promise and have opened up 
several lines of inquiry. In addition, the value of interactive modeling was shown, 
and a simple, fast model was developed which can be adapted in a matter of minutes 
to model any vibrating discrete systems. It could, if desired, be used to model wave 
motion in continua as well, constituting in that case a simple finite element model. 
All that is required is that the researcher correctly enter the equations of motion of 
the system to be studied, and to change the I/O section if needed to reflect those 
parameters only which are relevant to the new work. This moel may in the end 
prove to be the most important contribution of this research; in any case, it is hoped 
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that it will prove useful to many future researchers, to whom it is freely offered for 
modification and use. 
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APPENDIX A. NUMERICAL METHODS. 



Since all of the phenomena studied in this thesis have been oscillatory, it was 
possible to use a very simple numerical algorithm, which is a modification of the 
classical Euler’s method given by 

x n (t+h)-x n (t)+hv n (t), (107) 



and 

v n (t+h)-v n (t)+ha n (t). (108) 



where h is the time step used. Since the system which I intended to model was one 
where the rate of change of x depended at any given time on the values of x of the 
adjacent elements, it seemed reasonable to break the calculation of (1) up as follows: 

1. Calculate the acceleration felt by each element, based on the positions of 
its neighbors. 

2. Update the velocities of each element by multiplying the acceleration by 
the time step. 

3. After all velocities have been updated, calculate the new positions 
of the elements, and begin again at 1. 
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Mathematically, this process can be represented by a pair of equations: 




009) 



and 



x n (t+h)-x n (t)+hv n (t+h). 



( 110 ) 



This simple method turned out to work very well indeed. It was found that the 
free, linear lattice behaved as expected, and that the nonlinear driven lattice behaved 
as hoped. In order to check the quality of the numerics, an energy calculation was 
included in the program. The total energy of a free nonlinear lattice with linear 
coupling between nearest neighbors is given by 



It was found in all cases that the total energy of the system oscillated about some 
fixed value, with the amount of oscillation being directly proportional to the time 
increment used. Thus it was possible to use a large time step to observe qualitative 
changes in system behavior more expeditiously, and then to switch to a small step 







( 111 ) 



size (typically 0.005T, where T is the period of the oscillations) when precision is 



desired for quantitative comparisons to theory or other results. 
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As it turned out, this method and the conclusions about its efficacy had already 
been developed (Cromer [1981]), so I regretfully have to identify the method as the 
Euler-Cromer method... Also, the method solves a question posed by Goedde, 
Lichtenberg, and Lieberman [1990], who pointed out that, if too large a time 
discretization is used, parametric instabilities develop in the discrete Sine-Gordon 
Equation (of which our model is a third order approximation). They wondered, in 
their conclusions, whether it was possible to have stable solitons in a discrete Sine- 
Gordon system, or whether the energy of the system would always be 
equipartitioned among all of the elements when the time discretization was finite. 
The results presented in the body of this thesis, while not a rigorous proof, do 
strongly suggest that it is possible to have stable soliton solutions in a discrete 
system, provided the time step is reasonably small and a suitable method is used. 
Then, the errors are bounded, and one obtains stable results to at least many 
thousands of oscillator periods, so that even if the stability is only asymptotic, the 
model is practically stable. 

Whereas unfortunately the Euler-Cromer method was, if somewhat obscure, not 
original, the vehicle in which it was used is, I believe, original and of great potential 
utility to future researchers. In a sharp deviation from the type of numerical analysis 
which previous computer generations demanded, wherein one started the model with 
some set of initial conditions and parameter values and then waited for the outcome, 
repeating the process numerous times in order to build up a coherent picture of 
physical processes, this model system is highly interactive. This allows the researcher 
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to observe system behavior in "real time", and to adjust parameters as desired to 
control system behavior or to more fully understand it. Essentially, the model I have 
written is analogous to a complete laboratory setup, with data taking expedited and 
unwanted noise virtually eliminated (except numerical "thermal" noise, which is small 
in most cases). One can adjust the damping, drive amplitude, or drive frequency via 
either a fine or coarse adjust "knob", the knob being single keystrokes entered while 
the model is running at full speed . One may also "kick" any chosen elements by any 
desired amount, or "pin" an element at zero amplitude, in order to either perturb the 
system significantly or to move solitons along the lattice (as was done in the 
experimental studies with pendulum lattices and shallow water channels). 

The system’s behavior can be observed in the time domain either by watching 
the lattice on a full screen display of moving dots, or on a waterfall display where 
trends are more readily apparent (but detail is less clear). A zoom feature allows 
one to rescale with the lattice dynamics, or even to measure numerical "temperature" 
by zooming in on nodes in modes such as the " + 0-0" mode (see Chapter III), until 
the random noise caused by roundoff errors is detected. Several times this was 
done, and the results were a numerical noise so low that one had to zoom 2 40 times 
to see it! A user aid is available to monitor the system’s stability, which frees the 
user to some degree from having to continuously monitor the lattice. 

Additionally, the system’s dynamics can be viewed real time in the frequency 
and phase domains. In the former case, a running spectrum (consisting of 
consecutive 2048 point Fast Fourier Transforms adapted from Press, Teukolsky et. 
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al. [1988]) of a single element is displayed on the screen, along with all system 
parameters. In the latter, a phase portrait of up to five elements at a time is 
displayed, with Poincare section optional. In order to eliminate the drive frequency 
and thus clean up the display to make interpretation more straightforward, a rotating 
set of axes was used, with the coordinate and momentum axes rotating synchronous 
with the drive frequency. In such a display, a simple harmonic oscillator appears as 
a small ellipse (which vanishes as time increment approaches zero). 

The net result of the interactive features the model provides is that one is able 
to develop intuition about system behavior as various parameters are varied, and one 
can experiment with various types of perturbations to the system (including, as well 
as those just described, varying parameters randomly, via a gradient, or sinusoidally 
(Chapter IV), or just imposing a sudden ramp in amplitude of 5% to test stability of 
a given configuration). Such interaction leads one to find structures which might 
otherwise escape notice, since they often cannot be obtained simply by choosing a 
set of initial conditions and letting the system go. Additionally, having the variety 
of representations available at all times and in real time gives great diagnostic 
capability to the experienced researcher. 

Finally, the program was designed for maximum adaptability to the study of 
other problems concerning harmonic systems. In Appendix B, a detailed manual 
describes not only how to use the program, but how to alter it to model other 
systems. As discussed in Chapter IV, a working model of the Toda lattice was 
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available in fifteen minutes, with the complete range of interactive tools available 
for the study of that very challenging dynamic system. 
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APPENDIX B. PROGRAM MANUAL. 



A. INTRODUCTION AND GETTING STARTED. 

The program LATHCE.C is intended to be an adaptable modelling system for 
oscillating physical systems. The numerical method used, the Euler-Cromer method, 
is described in Appendix A, and should be understood by the researcher. The 
purpose of this appendix is to provide a readable manual which describes how to use 
and modify the program in order to study physical systems using the IBM PC 
(preferably 486 based, for speed). The program requires a VGA controller, and is 
intended to work in conjunction with screen-dump-to-printer programs such as 
EGALASER.COM, which is commercially available. The first part of this manual 
describes the use of the program, whereas the second is addressed to those who wish 
to make alterations to it in order to modify the model being studied. Changes simply 
to improve utility or to add functions are not addressed, although there are certainly 
many such changes which might be desirable. 

To start the program, simply type LATTICE at the DOS prompt. If you intend 
to use a data file for initial conditions, have its name, including extension, ready 
when you start the program. The program prompts the user for all information 
needed to start the program; interactive commands used once the program is started 
are not prompted. However, as discussed in the next section, the user may be 
prompted for additional information from time to time. If you desire to run the 
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program in conjunction with a screen dump program which is TSR (terminate and 
stay resident - it remains in your system’s RAM after you start it), you should start 
the TSR program before starting the model, since you cannot do it from within the 
model. 

B. USING THE PROGRAM. 

After you start the program, you will be asked whether you wish to use a data 
file for initial conditions. This is preferable, since otherwise you have little freedom 
concerning initial conditions (unless you modify the program yourself). If you do 
want to use a data file, type "y" or "Y" and < ENTER >. The program then asks 
whether you want to use the old format data file or not. The old format, which have 
been conventionally given the extension ".LAT", consist of a set of six parameters, 
followed by two numbers for each element. The first number for each element is 
its amplitude; the second is its velocity (which is usually close to zero). If this is the 
type of data file you have, type "y". Otherwise, your file should be formatted with 
only five initial parameters, but with each element having four numbers. These 
numbers represent, in order, the coupling, natural frequency, amplitude, and velocity 
of the element. Typing any character other than "y" at the prompt will result in the 
program trying to read this format. Finally, the program asks for the file name. The 
maximum size is 30 characters, so be sure the file name, including path if that is 
different from the path your program resides at, does not exceed thirty characters! 
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If you choose not to use a file for initial conditions, the program will prompt 
you for all of the parameters it needs. The program’s default initial conditions are 
an amplitude modulated cutoff mode. If the sign of the nonlinear coefficient is 
negative, corresponding to a hardening lattice, the program sets up an upper cutoff 
lattice; otherwise it sets up a lower cutoff lattice. The mode is modulated by one 
complete wavelength of amplitude chosen by the user ("modulation amplitude", as 
opposed to "mode amplitude", which is the amplitude of the cutoff mode itself). 

After either the file data is read or the model parameters have been obtained 
form the user, the program asks you to enter the time step. The number you enter 
will be the time step, after it is multiplied by one half per cent of the period of the 
drive (this requires that you do not enter zero for drive frequency , even if you want 
an undriven system, since doing so will give a divide by zero error, and this program 
doe, not include error checking !. Typical choices for time step, as entered, will be 
from one half to four, depending on whether you are most interested in exact 
amplitudes (use small time step) or in getting relatively quick, qualitatively accurate 
results (use large time step). Use of time steps greater than four does not cause an 
error, but it should be done only after the user is familiar with the behavior of the 
model he is studying. If the model is near any sort of boundary in behavior type, or 
if it is near a potential energy maximum, using a large time step may perturb the 
system into undesirable or even catastrophic behavior. On the other hand, judicious 
use of very large time steps can push the system into stable behaviors which might 
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otherwise have been missed — experience is important when choosing the time step! 



1. The Default Graphics Mode. 

At this point, the screen will go blank. Type ~G to get a graphics display 
of the system in real time (^G means < CONTROL > G, and this type of shorthand 
will be used from now on). The ~G command will always do this, so it is a useful 
one to memorize immediately. If the display is going off the screen due to excessive 
amplitude, you can "zoom" out by typing "o". This means in effect that you step 
back from the system, so that the observed maximum amplitude increases (by a 
factor of two, as it happens). Similarly, "i" zooms you back in by a factor of two. 
This is useful if the amplitudes get too small to see (your picture is limited by the 
pixel resolution of the graphics mode). The graphics display will only show 40 
elements on the screen; if the lattice you are using contains more than that, the forty 
visible elements are the middle elements. The other elements are still being 
calculated and will be included if you save your state to a data file, you simply can’t 
see them. The maximum number of elements is 150. 

As a user aid, a stability checking routine is included in the program. This 
routine, which is called by typing "S" from any point in the program, asks the user 
to select an element for monitoring and then monitors that element to see if it is 
"stable". When choosing an element, remember that, in the C language, the first 
element of an array is the zeroth element. If you want to monitor the first element, 
tell the computer to monitor element zero. Stability is determined by comparing 
successive peaks of the amplitude of the chosen element. If twenty peaks in a row 
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are all within one percent of each other, the system is called stable. This does not 
mean it is stable in any scientific or mathematical sense; this function is solely to aid 
the user as he sees fit. When this option is chosen, the entire lattice changes color 
to red, except the element being monitored (which remains white). After stability 
has been verified, the lattice returns to the white color, except the chosen element, 
which then becomes red. 

The stability checker is one of two routines that uses the amplitude peak 
detection routine -- the other is the Waterfall Display, which is described later. The 
stability checker must be active if it is desired to save the current model state to a 
data file, since the states are saved only when the monitored element reaches an 
amplitude peak. If, while the stability checker is active, you type "p", the program 
will pause and ask you whether you wish to save the state to a file. If you do, type 
"y" or "Y", and then give the file name. It is recommended that a personal file 
naming protocol be developed early so that confusion does not occur when the 
program suddenly asks for a file name. If you type "p" before the stability checker 
is inactive, the pause routine will be executed as soon as the stability checker is next 
activated. After dumping data to a file, or declining to do so, the program prompts 
for a new time step. The old one is displayed for convenience. If you want to keep 
it, you must type it. The program uses whatever number you type as the new time 
step. Thus, the pause routine provides the only means for changing the time step 
during model operation. The user is cautioned against changing the time step while 
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the Fourier Analysis routine is operating (see below), since a change in the time step 
then will corrupt the Fourier data and give erroneous spectra. 

The drive parameters and the dissipation parameter can be changed while 
the program is operating (it does not need to be in Default Graphics Mode to do 
this). The letters U,V,D, and E represent coarse and fine increases and decreases 
in the applicable parameter. A <SHIFT> operation with one of these letters means 
the parameter affected is drive frequency. Thus U means coarse increase in drive 
frequency. Similarly, means coarse increase in drive amplitude, and u means 
coarse increase in dissipation. The commands ^Z,Z, and z mean that the indicated 
parameter (drive amplitude, frequency, and dissipation, respectively) gets set directly 
to zero. This is very useful if your model starts to blow up -- type ~Z and eliminate 
drive, so that system energy will at best remain constant (it will decrease if there is 
any damping still present). Finally, s brings all elements immediately to rest at zero 
amplitude (another useful, but destructive, way to stop a system from blowing up). 

It is possible to interact on an element by element basis with your model 
system, using the and keys. The first of these "kicks" an element, which 
means a step change in amplitude is made. Velocity does not change, however, 
which should be kept in mind by the observer. The second "pins" an element at zero 
amplitude and velocity. The element stays pinned until is typed again and the 
same element selected (i.e., is a toggle switch, as are several other commands). 

The lattice can be changed on a global scale by introducing 
nonuniformities into the parameters representing coupling and natural frequency. 
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using the W keystroke. These parameters are stored as arrays, with a separate value 
for each element. This is so they may be made to vary from element to element. 
There are three types of variation available, and they can be applied successively as 
many times as desired (this includes "undoing" a variation by introducing an identical 
but opposite polarity variation -- except, of course, for random variations). The first 
type is a random, five percent (peak to peak) variation of the parameter (only one 
of the two parameters can be varied for each keystroke, although both can be varied 
by successive keystrokes). The second type is a linear gradient. This gradient is 
positive, with the zeroth element’s parameter remaining unchanged, and the last 
element’s parameter being changed by the full amount of the gradient. This amount 
(in percent) is chosen by the user. The third type is a sinusoidal variation, in which 
the amplitude of the sinusoidal variation and the (integer) number of wavelengths 
of the modulation along the entire length of the lattice are input by the user. These 
parameter variations are useful for studying the effects of nonuniformities on the 
dynamics of the lattice being studied. Additional types of variations can easily be 
added by the user (see modification instructions below). 

Another useful function, dumping a single element’s behavior to a data file 
in the form of a pair of time series representing displacement and velocity, is 
available to the user in all modes ("modes" refers to Default Graphics, Phase Plot, 
Text, Waterfall, and Spectrum display modes). Type ~F to start file dumping. The 
program will prompt for the element to be dumped; then the program returns to 
normal operation. To stop the dump, type *F again, and give the name of the data 
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file you want to output to. If you don’t type again, the program will prompt you 
after the 4096 point buffer is full. 

Finally, typing n perturbs the system by kicking each element by up to 
10%, depending on its location. This command is useful when the stability to 
perturbation of a state is to be determined. If the state returns to normal after a 
few of these perturbations, it is definitely stable! 

2. Phase Plot Mode. 

Phase plane diagrams, or phase plots, are very useful tools in gaining an 
understanding of dynamical systems. Accordingly, this type of display is available 
to the user via the P command. This command, which can be used anytime the 
model is running, shifts the display to phase plot mode. The program asks whether 
Poincare sections are desired. Typically, they are, since they provide a much cleaner 
and easier way to interpret display. Type y to get Poincare sections, anything else 
to get phase plots with continuous updates (one dot per time step). The program 
then asks for up to five elements to display. Entering five will result in the model 
resuming and the phase plot beginning to build up. Entering fewer than five results 
in nothing happening until the user enters 999 to indicate that no further elements 
are desired (this entry is prompted, so you don’t have to remember it). When in 
phase plot mode, all of the commands which affect the lattice are valid, and their 
effects will be seen in the phase plots. It is recommended that you experiment with 
this a lot, because it is very interesting and useful, and phase plot information is 
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file you want to output to. If you don’t type again, the program will prompt you 
after the 4096 point buffer is full. 

Finally, typing n perturbs the system by kicking each element by up to 
10%, depending on its location. This command is useful when the stability to 
perturbation of a state is to be determined. If the state returns to normal after a 
few of these perturbations, it is definitely stable! 

2. Phase Plot Mode. 

Phase plane diagrams, or phase plots, are very useful tools in gaining an 
understanding of dynamical systems. Accordingly, this type of display is available 
to the user via the P command. This command, which can be used anytime the 
model is running, shifts the display to phase plot mode. The program asks whether 
Poincare sections are desired. Typically, they are, since they provide a much cleaner 
and easier way to interpret display. Type y to get Poincare sections, anything else 
to get phase plots with continuous updates (one dot per time step). The program 
then asks for up to five elements to display. Entering five will result in the model 
resuming and the phase plot beginning to build up. Entering fewer than five results 
in nothing happening until the user enters 999 to indicate that no further elements 
are desired (this entry is prompted, so you don’t have to remember it). When in 
phase plot mode, all of the commands which affect the lattice are valid, and their 
effects will be seen in the phase plots. It is recommended that you experiment with 
this a lot, because it is very interesting and useful, and phase plot information is 
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most useful only after significant experience has been built up looking at them under 
various conditions. 

3. Text Mode. 

The text mode is activated by 'T. If this is preceded at any time by ~H, 
the text mode display will conclude with a display of total system energy', calculated 
as described in Appendix A. In this mode, the text mode functions as a mode in 
which model operation continues; in any other case, the model operation is 
suspended until you leave text mode by using the appropriate command to enter 
another mode. The text mode gives the values of all system parameters, and then 
waits for you to hit a key. When you do, it displays the amplitude and velocity of 
each element. In all honesty, this is not a very refined mode, so it is really useful 
only for looking at system parameters and system energy. 

4. Spectrum Mode. 

By typing ''X, the user enters Spectrum Mode. The user is prompted to 
enter the number of the element whose spectrum is desired. Then execution returns 
to the mode the user was in when he typed ^X. However, the amplitude of the 
chosen element is now dumped every tenth time increment to a data array of 2048 
elements. When the array is full, the Fast Fourier Transform is calculated and 
displayed. This display then remains on the screen until it is updated (the model 
continues to run, and sequential FFTs are calculated every 20480 time increments 
until the user enters a different mode using the appropriate command). The system 
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parameters and maximum frequency shown are displayed at the top of the screen, 
and a frequency axis is displayed at the bottom so that numerical values of the 
spectral components can be estimated. If it is desired for the complex FFT values 
to be dumped to a data file for later analysis, type y and then enter the filename 
when prompted. 

5. Waterfall Mode. 

The waterfall mode, entered via the ~W command, is useful for 
monitoring trends, particularly after some sort of change has been made. I found 
it most useful when I introduced lattice nonuniformities. If I immediately typed ~W, 
I got a good visual indication of the exact response of the lattice to the step changes 
I introduced. This type of use is commended to the user. When in Waterfall Mode, 
a waterfall display of 1 10 iterations is displayed, one iteration at a time. It is best 
to enter Waterfall Mode with the stability checker inactive, to avoid multiple 
waterfall entries for the same lattice state. If this is done, the lattice is displayed 
once during each cycle, at the same point in the cycle, so that useful comparisons can 
be made. When the Waterfall Display is full, the system does not continue to 
operate . This is so that one can pause to examine the display, dump it to the printer 
if desired, and then press ^W again to restart the model and refresh the display. 
Operating in this way makes it possible to have a continuous series of waterfall 
displays recording system behavior for long periods of time. 

The zoom in and out commands (i and o) are operative in this mode (and 
the phase plot mode), so the display can be adjusted for maximum effectiveness. 
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Also, typing ~W while the waterfall display is not full toggles the display off. The 
model then continues to run, but no new display will appear until the appropriate 
mode command (''G, ~T, P, ~X, or *W) is entered. 

When finished with the program, type ~Q to quit. 

C. MODIFYING THE PROGRAM. 

The program was designed to be as modular as possible, although there was no 
great effort expended to optimize performance in general, as the emphasis was on 
research and not programming finesse. This modularity makes it easy to modify the 
program to model other physical systems. Other modifications, such as to add 
additional features or display modes, will not be addressed here. 

To model a physical system, the first requirement is to develop the discrete 
equations of motion. Discreteness is required for computational purposes, although 
time can be treated as continuous (it is discrete, but the equations do not need to 
reflect that, since we calculate the acceleration and velocity at discrete instants in 
time, thus implicitly discretizing time for you). These equations should then be 
coded up in the C language, and then inserted directly in a copy of the waterfall 
program in place of the equations of motion currently installed. The equation 
entered should be of the form 
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x -^coordinates, momenta,, f) . 



( 112 ) 



As long as this form is used, and as long as the system described by the equations 
of motion is oscillatory, at least in steady state, then the numerical methods used 
should provide excellent results. The user should verify this in every case by 
checking the numerics as described in Appendix A. 

Once the equations of motion are modified, you need to initiate any variables 
or parameters that appear in the equations of motion that are not already initiated. 
For example, when I modified the program to create a Toda Lattice model (see 
Chapter IV), I had to add a variable b, representing the exponential coupling 
parameter. If these parameters are to be input by the user, then the user_init 
function should be modified; also, any references to parameters which are removed 
in the new model should be deleted from the I/O statements to minimize confusion. 
The parameter changes should also be reflected in the file I/O sequences in 
user_init and process_pause functions. Also, for completeness, the user should 
modify the parameters displayed in display_text, init_phasplot, and plot spectrum 
functions. When these changes are made, the model is ready to compile using either 
the QUICKC compiler or, if time optimization is required, the Microsoft C version 
6.0 compiler. That is all there is to modelling a new physical system! As I stated 
in Chapter IV, this is a real strength of this program; I created a Toda lattice model 
in fifteen minutes. 
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SUMMARY OF COMMANDS 



~Q 

"G 

Arjp 

"W 

P 

"X 



p 

w 

~H 

"F 

i/o 

"L 

"K 



y 

n 

S 

s 

"on 

"U,"V,"D,"E 

U,V,D,E 

u,v,d,e 



Quit the Program. 

Enter Default Graphics Mode. 

Enter Text Mode. 

Enter Waterfall Mode. 

Enter Phase Plot Mode. 

Enter Spectrum Mode. 

Pause program (file I/O, change time step). 

Introduce nonuniformities into lattice. 

Monitor total system energy (in text mode only). 

Toggle time series data dump to file. 

Zoom in/out by a factor of two in amplitude. 

Pin an element. 

Kick an element. 

Dump spectrum to data file. 

Pertum the system by a ten percent amplitude gradient. 
Toggle stability checker on/off. 

Stop all elements. 

Increase/decrease Poincare strobe frequency. 
Coarse/fine increase/decrease of drive amplitude. 
Coarse/fine increase/decrease of drive frequency. 
Coarse/fine increase/decrease of dissipation. 
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~Z,Z,z 

~R 



Zero out the drive amplitude/frequency/dissipation. 
Reset and restart the program. 
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C. PROGRAM CODE 



/* PROGRAM WATERFALL GENERALIZED NONLINEAR LATTICE MODEL 

VERSION 3.0 

WRITTEN BY BRIAN GALVIN 
LAST_UPDATE 10 JAN 1991 

THIS PROGRAM SIMULATES A GENERAL LATTICE WITH EQUATIONS OF MOTION 
THAT CAN BE SUBSTITUTED IN WHERE INDICATED. VARIOUS INTERACTIVE 
FEATURES ARE PROVIDED , WHICH ARE EXPLAINED IN THE PROGRAMMER ' S 
MANUAL. A WATERFALL DISPLAY IS ADDED TO MONITOR SOLITON MOTION 
DUE TO MEDIUM NONUNIFORMITY EFFECTS. COUPLING AND NATURAL 
FREQUENCY CAN BE MADE NONUNIFORM VIA SHIFT-W. */ 



/* *** PREPROCESSOR DIRECTIVES *** */ 



linclude 

linclude 

linclude 

linclude 

linclude 

linclude 

linclude 



"\quickc\inc\BIOS .H M 
M \quickc\inc\graph.h" 
M \quickc\inc\stdl ib.h" 
"\quickc\inc\CONIO. H M 
"\quickc\inc\MATH . H M 
M \quickc\inc\STDIO . H M 
M \quickc\inc\time.h M 



/* *** MACRO DEFINITIONS *** */ 

Idefine SQR(a) ((a)*(a)) 

Idefine CUB(a) ( ( a ) * ( a ) * ( a ) ) 

Idefine SWAP(a,b) tempr=( a) ; (a )=(b) ; (b)=tenpr /* USED IN FFT ROUTINE */ 

Idefine DOFOR(i,to) f or ( i=0 ; i <to ; i ++ ) 

Idefine PI 3.14159265359 
Idefine SCREEN_CORRECTION_FACTOR 1.05 
Idefine eta_increment 0.01 
Idefine beta_increment 0.01 
Idefine oroega__increroent 0.001 
Idefine STABILITY_INCREMENT 0.01 
Idefine DISPLAY_INCREMENT DINC 
Idefine PHASE_INCREMENT PINC 
Idefine MIDDLE_ELEMENT ( no_pendulums/2 ) 

Idefine text_flag flags[0] /* flags have names in body of program, */ 

Idefine graphics_f lag flags[l] /* but they all form a single array. */ 

Idefine f ile_dump_f lag flags[2] 

Idefine stability_f lag flags[3] 

Idefine peak_flag flags[4] 

Idefine stable_flag flags[5] 

Idefine phase_flag flags[6] 

Idefine pause_flag flags(7] 

Idefine pause_flag2 flags[8] 

Idefine stop_flag flags[9] 

Idefine spectrum_f lag flags[10] 

Idefine dump_spectrum_f lag flags(ll) 

Idefine energy_flag flags[12] 

Idefine waterf a 1 l_f lag flags[13] 

Idefine wait_flag flags[14] 

/* *** DYNAMICAL VARIABLE DECLARATION. THESE ARE GLOBAL VARIABLES *** */ 

double coordinate [ 150 ] , momentum [ 150 ] , old_coordinate [ 1 50 ] ,old_momentum[ 150 ] ; 



/* SIMPLIFIED DO LOOP COMMAND */ 

/* These increments are used in */ 

/* adjusting parameters */ 



double acceleration , gamma [ 150 ] , roean_g aroma , beta , omega 0 [ 150 ] ,mean_oroegaO , 
eta , omega , alpha , model _ti me , time_int , 

max_amp , mode_amp , t iroe_series_disp[ 8000 ] , phase_elements[ 5 ] , 
peak_record[ 20 ] , pinned_elements[ 150] , DINC, PINC, period. 
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spectrum [ 4098 ] ,temp2[4096 ] , energy , gradient ,wat_inc; 
int no_pendulums , f lags[ 15 ] ,counterl ,phase_counter ,phase_int , f irst_eleroent , 
chosen_element ,stability_eleroent ,disp_co lor , colors [ 5 ] , counter 2 , 
spectrum_element , spectrum_counter , sample_counter , energy_counter , 
waterf all_counter , color_counter ; 



main() ( 

FILE *fp; 

int c , i , j ,k, 1 ,m,n ,xx; 
double modulation; 
char ans [ 5 ] ; 

_clearscreen(_GCLEARSCREEN ) ; 

user_init ( ) ; /* STARTS PROGRAM WITH USER INPUT */ 

_setvideomode(_MRES4C0L0R) ? 

_setcolor ( 1 ) ; 
xx«0 ; 

wat_inc=l ; 
disp_color=l ? 

/* INITIALIZE COUNTERS HERE */ 

counter 2=0; 
energy_counter=0 ; 
color_counter=o ; 



/* *** BODY OF PROGRAM STARTS HERE *** */ 



while ( stop_flag==0 ) ( 
energy=0; 

if (pause_f lag2==l ) { 
pause_f lag2=0 ; 
pause_f lag=0 ; 
process_pause( ) ; 

} /* if ( pause . . . ) */ 

if (kbhit( ) !«=0) ( 
c=getch( ) ; 
if ( C “17 ) ( 

stop_f lng~l ; 

) 

else if (c==112) ( 

pause_f lag=l ; 

) 



/* A peak has been detected, so the */ 

/* pause routine is called */ 

/* CHECK FOR USER INPUT DURING MODEL RUN */ 
/* *Q : Quit program */ 

/* p : Pause program/dump state */ 



else if ( c~20 ) ( /* *T : Text Mode ♦/ 

_clearscreen(_GCLEARSCREEN ) ; 
text_f lag=l ; 
graphics_f lag=0 ; 
waterf all_f lag=o ; 
wait_f lag=0 ; 
spectrum_f lag=0 ; 
display^ text( ) ; 



) 

else if(c==7) ( /* A G : Graphics Mode */ 

_clearscreen(_GCLEARSCREEN) ; 
_setvideomode(_MRES4COLOR) ; 

_setcolor ( 1 ) ; 
phase_f lag=0 ; 
wait_f lag=0; 
spectrum_f lag=0 ; 
watorfal l_f lng-0; 
text_f lag=0 ; 
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graphics_f lag=l ; 

) 

else if(c==8) ( /* ~H : Monitor energy in text mode */ 

if (energy_f lag==0) energy_f lag=l ; 
else energy_f lag=0 ? 



/* *** THE FOLLOWING CASE IS FOR CHANGING THE UNIFORMITY OF THE 
LATTICE. COUPLING OR NATURAL FREQUENCY CAN BE VARIED. 

THE THREE TYPES OF VARIATION ARE: 

1. RANDOM VARIATION OF FIVE PER CENT 

2. LINEAR GRADIENT OF ANY AMOUNT 

3. SINUSOIDAL GRADIENT OF ANY AMOUNT, WITH AN 

ARBITRARY NUMBER OF FULL WAVELENGTHS OVER 
THE LENGTH OF THE LATTICE *** */ 

else if(c=~87) ( /* W : Kick in parameter variations */ 

_setvideomode(_DEFAULTMODE) ; 
srand( (unsigned) time(NULL) ) ; 

printf ( "\nEnter 1 if you wish to vary coupling: " ) ; 

scanf ("%d",& j); 
if(j~l) ( 

printf ( "\nEnter 1 (random), 2 (gradient), or 3(sine) desired: ") 

scanf ( "td M , fck ) ? 
if ( k==2 ) ( 

printf ( "\nEnter gradient (percentage over entire lattice): H ); 

scanf ( "%lf " , ^gradient ) ? 

DOFOR( i ,no_pendulums ) ( 

gamma( i )=gamma( i ]+mean_gamma*i*gradient/( 100*no_pendulums ) : 

) /* DOFOR */ 

) /* if ( j“2 ) */ 
else if ( k==3 ) ( 

print f ( "\nEnter modulation amplitude (percentage): " ) : 
scanf ( n \ 1 f *' , ^gradient ) ; 

printf ( "\nF.nter integer number of wavelengths to use: " ) ; 

scanf ( "Id" , & 1 ) ; 

DOFOR ( i , no_pcndu 1 urns ) ( 

gamma [ i ]=gamma( i ) +mean_gamma* gradient/1 00* ( -cos( 2*l*PI*i/ 
no_pendulums ) ) ; 

) 

) 

else DOFOR( i ,no_pendulums) < 
j=rand( ) ; 

gamma ( i J^mean^gamma* ( . 95+ . 1 * j/RAND_MAX ) ? 

) 

) 

else ( 

printf ( M \n Varying omegao. . . \n" ) ; 

printf ( "\nEnter 1 (random), 2 (gradient), or 3(sine) desired: " ) 

scanf ( "%d M , fik) ; 
if ( k==2 ; ( 

printf ( *'\nEnter gradient (percentage over entire lattice): "); 

scanf ( •• % 1 f •• , ^gradient) ; 

DOFOR( i , no_pendulums ) ( 

omegaOf i ]=omega0[ i J+mean^omegoO* i *gradient/( 1 00*no_pendulums ) ? 

) 

) 

else i f ( k==3 ) ( 

printf ( "\nEnter modulation amplitude (percentage): " ) ; 

scanf ( "%lf ” , ^gradient ) ; 
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printf ( "\nEnter integer number of wavelengths to use: "); 

scanf (”%d",61) ; 

DOFOR ( i , no_pendu 1 urns ) { 

omegaO[ i ]=oraegaO[ i ]+mean_omega0*gradient/100* ( -cos ( 2*l*PI*i/no_pendulums ) ) ; 

> 

) 

else ( 

DOFOR ( i , no_pendulums ) { 
j=rand ( ) ? 

printf ( "Random number is: %d RAND_MAX is %d\n M , j ,RAND_MAX) ; 
omegaO [ i ] =mean_omegaO* ( . 9 5+ . 1 * j/RAND_MAX ) ; 

) 

) 

) 

) 

else if(c==23) { /* A W : Waterfall Display Mode */ 

if (waterf al l_f lag==0 ) { 
phase_f lag=0 ; 
graphics_f lag=o ; 
text_f lag=0? 
spectrum_f lag=0 ; 
waterf a 1 l_f lag=l ; 
init_waterfall( ) ; 

) 

else { 

if (wait_f lag==l ) ( 
wait_f lag=0 ; 
init_waterf all( ) ; 

) 

else waterf al l_flag=0 ; 

) 

) 

else if(c==80) ( /* P : Phase Plot Mode */ 

if (phase_f lag==o) init_phasplot ( ) ? 
else { 

phase_f lag=0 ; 
spectrum_f lag=0 ; 
graph ics_flag=l ; 
waterf all_f lag=0 ; 
wait_f lag=0 ? 

_setvideomode(_MRES4COLOR) ; 

_clearscreen(_GCLEARSCREEN) ; 

) /* else */ 

) 

else if(c==6) { /* *F : Time series to file */ 

if ( f ile_dump_f lag==0 ) start_f ile_dump( ) ; 
else { 

stop_f ile_dump( ) ; 
f ile_dump_f lag=0 ; 

) 

) 

else if ( c==24 ) ( /* *X : Calculate spectrum */ 

spectrum_f lag=l ? 

_setvideomode(_DEFAULTMODE) ; 

_clearscreen(_GCLEARSCREEN) ; 

printf ( M \nEnter number of element to be analyzed: " ) ; 

scanf ( "%d" ,&spectrum_element ) ; 

_setvideomode(_MRES4COLOR) ; 

_clearscreen(_GCLEARSCREEN) ; 
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else if (c==21 ) ( * /* *U : 

eta=eta+eta_increment; 
stability_restart ( ) ; 

) 

else if(c»»4) ( /* *D : 

e t a**e t a -eta_ increment; 
stability_restart( ) ; 

) 

else if (c er =22 ) ( /* A V : 

eta®eta+eta_increment/10 ; 
stability_restart ( ) ; 

) 

else if(c==5) ( /* *E : 

eta=eta-eta_increment/10 ; 
stability_restart ( ) ; 

) 

else i f ( c==26 ) ( /* *S : 

eta=0; 

stability_restart( ) ? 

) 

else if (0**85 ) ( /* U : 

omega=omega+omega_increment ; 
stability_restart( ) ; 

> 

else if(c=-68) ( /* D : 

time_int-tiroe_int* 200/period; 
omega=omega-omega_increraent ; 
period=2*PI/omega ; 
time_int=time_int*period/200 ; 
stability_restart( ) ; 

) 

else i f ( c==86 ) ( /* V : 

time_int=time_int* 2 00/period; 
omega=omega+omega_increment/l 0 ; 
period=2*PI/omega ; 
time_int=time_int* period/ 200; 
stabi lity_restart ( } ; 

) 

else if ( c==69 ) ( /* E : 

time_int=time_int* 200/period ; 
omega =-omcga-omega_increment; 
per iod=2* PI /omega ; 
time_int=time_int*period/200 ; 
stability_restart( ) ; 

) 

else if ( c==90 ) ( /* s 

omega =0 ; 

stability_restart( ) ; 

) 

else if (c==117) { /* u 

beta=beta+beta_increment ; 
stability_restart ( ) ; 

) else if(c= 

beta=beta-beta_increment ; 
stability_restart( ) ; 

) 

else if (c==l 18 ) { /* v 

beta=beta+beta_increroent/10 ; 
stability_restart ( ) ; 

) 



Increase drive (coarse) */ 



Decrease drive (coarse) */ 



Increase drive (fine) */ 



Decrease drive (fine) */ 



Turn off drive */ 



Increase frequency (coarse)*/ 



Decrease frequency (coarse)*/ 



Increase frequency (fine)*/ 



Decrease frequency (fine)*/ 



Set frequency to zero */ 



Increase damping (coarse)*/ 

=100) ( /* d : Deer 



: Increase damping (fine)*/ 
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else if(c==121) ( /* y : Dump spectrum */ 

if (spectrum_f lag==l ) dump_spectrum_f lag=l ; 



) 

else if ( c~=101 ) { /* e 

beta«beta-beta_increroent ; 
stability_restart( ) ; 

) 

else if ( c“122 ) ( /* z 

beta=0 ; 

stability_restart( ) ; 

) 

else if ( c*=12 ) ( /* 

pin_elements ( ) ; 
stability_restart( ) ? 

) 

else if (c==ll) { . /* *K 

kick_elements( ) ; 
stability_restart( ) ; 

) 

else if ( c==105 ) { /* i 

DINODINC/2; 

PINOPINC*2; 

wat inc=wat_inc*2 ; 

) 

else if (c==lll ) ( /* o 

DINC=DINC*2; 

PINC=PINC/2? 
wat_inc=wat_inc/2 ; 

) 

else if (c==110) ( /* n 

DOFOR( i ,no_pendulums) ( 

modulations l*i/no_pendulum 
coordinate[ i ]=coordinate[ i ) 



: * Decrease damping (fine)*/ 



: Turn off damping */ 



: Pin elements */ 



: Kick elements */ 



: Zoom in */ 



: Zoom out */ 



: Perturb system */ 

( 1 . OO+modulation) ; 



) 



) 

else if ( c==15 ) ( /* *0 : Decrease strobe freq */ 

phase_int=phase_int-l ; 

) 

else if(c==9) ( /* *1 : Increase strobe freq */ 

phase_int=phase_int+l ; 

) 

else if(c==18) ( /* *R : Restart program */ 

_setvideomode(_DEFAULTHODE) ; 

_clearscreen(_GCLEARSCREEN) ; 
user_init( ) ; 
stability_restart ( ) ; 

if ( phase_f lag==0 ) _setvideomode(_MRES4COLOR ) ; 
else _setvideomode(_VRES16C0L0R) ; 



) 

else if ( c==115 ) ( /* s 

DOFOR( i ,no_pendulums) ( 
coord inate[ i ]=0; 
momentum [ i ]=0? 

) /* DOFOR */ 
stability_restart( )? 

) 

else if(c==83) ( /* S 

monitor_stability ( ); 

) 

/* if ( kbhit . . . ) */ 



Stop everything as is */ 



Monitor stability */ 
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spectrum_counter++ ; 
sample_counter=0 ; 

} 

else | 

compute_spectrum( ) ; /* COMPUTE SPECTRUM IF SPECTRUM DATA FULL ♦/ 

_clear screen ( _GCLEARSCREEN ) ; 

_setvideoraode(_VRES16COLOR) ; 
plot_spectrum( ) ; 
phase_f lag=0; 
text_f lag=0; 
graphics_f lag=0 ; 

DOFOR (j,4096) spectrum[ j ] =0 ; 
spectrum_counter=0 ; 
sample_counter=0 ; 

} /* else */ 

} /* if ( (spec. . . */ 

/* CHECK FOR PEAKS AND MONITOR STABILITY/SET PAUSE FLAG/INITIATE 
WATERFAL DISPLAY AS REQUESTED BY THE USER ♦/ 

if ( ( (i=^stability_element)&6(stability_f lag==l ) && ( stable_f lag==0) ) | 

( ( i==(no_pendulums-l ) )&& (waterfall_f lag==l ) ) ) { 

if ( ( coordinate [ i ]>0 )&& ( coordinate ( i ]<old_coordinate( i ] ) 

&& ( peak_f lag==0 ) ) ( 
peak_f lag=l ; 

if (pause_f lag==l ) pause_f lag2=l ; 

DOFOR (k, 19 ) peak_record[k ]=peak_record[ k+1 ] ; 
peak_record( 19 ]=coordinate[ i J ; 
if ( stability_f lag==l ) stabili ty_check( ) ; 
if ( waterf all_f lag==l ) disp_waterf all ( ) ; 

) 

else if (coordinated J<0) peak_flag=0; 

) 

) /* DOFOR */ 

®odel_time=inodel_tiine+tiine_int ; /♦ INCREMENT TIME */ 

/* CHOOSE DISPLAY MODE AND EXECUTE IT ♦/ 

if (graphics_f lag==l ) display_graphics( ) ; 

if ( (text_f lag==l )&&(energy_flag==l ) ) ( 

printf("\n Total Energy at time %lf is %lf” , model _ti me , energy ) ; 
energy_counter++ ; 

^ if (energy_counter==20 ) ( 

printf ( "\n\nStrike any key to continue...”); 
while(kbhit( )~0) ; 
energy_counter=0 ; 

) 

) 

if ( phase_f lag !-0) ( 
phasplot( ) ? 

phase counter++; /* Used for Poincare sections */ 

} 

) /* if(wait_flag ♦/ 

) /* while(stop_f lag. . . */ 

/* ««* LEAVE MAIN BODY OF PROGRAM HERE ♦/ 
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/♦ *** EQUATIONS OF MOTION. THE FOLLOWING LOOP CALCULATES THE 
ACCELERATION OF EACH ELEMENT AS THE FIRST STEP OF THE 
EULER-CROMER METHOD. THE EQUATIONS CAN BE CHANGED TO 
MODEL DIFFERENT PHYSICAL SYSTEMS WITH NO OTHER CHANGE 
NEEDED (ALTHOUGH IN ALL LIKELIHODD THE USER_INIT AND 
FILE SAVE ROUTINES SHOULD BE ALTERED IF THE PARAMETERS 
OF INTEREST HAVE CHANGED. ♦♦♦ ♦ /. 

if (wait_f lag”0) ( 

DOFOR ( i , no_ pendulums ) ( ’ 
if(i=0) ( 

acceleration=gamma [ i ]♦ (coordinate! i + 1 ]+coordinate [no_pendulums-l ] 
-2*coordinate[ i ] ) -bet a ♦momentum! i ] 

-( SQR(omegaO[ i ] ) +2*eta*cos( 2 ♦omega ♦modest ime ) ) ♦coordinate! i ] 
+alpha*CUB(coordinate[ i ] ) ? 

) 

else if ( i“( no_pendulums-l ) ) ( 

acceleration=gamma[ i ]♦ (coordinate! i-1 ]+coordinate[ 0 ] 

-2 ♦coordinate! i ) ) -beta ♦momentum [ i ] 

-( SQR ( omega 0 [ i ] ) +2*eta*cos( 2*omega*model_time ) ) ♦coordinate! i ] 
+alpha*CUB( coordinate! i ] ) ; 

) 

else ( 

acceleration=gamma [ i ]♦ ( coordinate! i+1 ]+coordinate ( i-1 ] 

-2 ♦coordinate [ i ] ) -beta ♦momentum! i ] 

- ( SQR( omegaO l i ] ) +2*eta*cos( 2*omega*model_time ) ) ♦coordinate! i ] 
+alpha*CUB(coordinate[ i ] ) ; 

) 

old_momentum! i ]=momentum[ i ] ; 

momentum! i ]=momcntum[ i ] +acccl erat ion ♦time. int ? 

) /♦ DOFOR ♦/ 

/♦ UPDATE ACTUAL ELEMENT POSITIONS BASED ON THE NEW VELOCITIES ♦/ 

DOFOR ( i ,no_pendulums ) ( 

old_coordinate[ i ^coordinate ( i ] ; 
if (pinned_elements[ i ]==0 ) ( 

coordinate! i ^coordinate [ i ]+monentum! i ]*time int; 

) 

if ( ( i~=chosen_element ) ( f ile_dump_f lag==l ) ) ( /♦ DUMP TO FILE BUF */ 
if (counter2++==30 ) ( 

time_series_disp[counterl++ ]=coordinate[ i ] ; 
counter2=0 ; 

) 

if (counterl==8000) stop_f ile_dump( ) ; 

) /♦ if ( ( i — ♦/ 

if (energy_f lag~l ) ( /♦ CALCULATE ENERGY IF REQUESTED ♦/ 

energy=energy+0 . 5* ( SQR( momentum! i ] ) 

+ gamma[ i ]*SQR( (coordinate! i+1 ] -coordinate! i ) ) ) 

+SQR( coordinate! i ] ) ) -a lpha/4^pow( coordinate! i ] , A ) ; 

) /♦ if (energy ... ) ♦/ 

if ( (spectrum_flag"l)&&(spectrum_element«==i) /♦ UPDATE SPECTRUM DATA */ 
&& ( (sample_counter++ )=-5 ) ) ( 
if ( spectrum_counter<2048 ) ( 

spectrum! 2*spectrum_counter Incoordinate! i ] ; 
spectrum! 2+spectrum_counter+l 1=0; 



196 



_setvideomode(_DEFAULTMODE) ; 

printf ("PROGRAM COMPLETE AT \lf " , model_time ) ; 

) /* MAIN ( ) */ 

/***♦ USER_INIT ACCEPTS USER INPUT FOR STARTING PARAMETERS AND INITIAL 
CONDITIONS , INCLUDING CHOICE OF DATA FILE TO START FROM. ALSO, 
VARIOUS INITIALIZATIONS ARE PERFORMED TO GET THE SYSTEM GOING **♦*/ 

user_init() ( 

FILE *fq; 

int c,i,j,k,l,m; 

char answer[ 1 ] , f ilename[ 30 ] ; 

printf ("GENERALIZED LATTICE MODEL PROGRAM W/ VARIABLE PARAMETERS^"); 

printf ("Version 2.1X486 (MS) \n"); 

printf ("Last updated 10 JAN 1991\n" ); 

printf ( "Variant notes: Default IC is AM\n" ) ; 

printf ( "OmegaO is set equal to one for all cases !\n"); 

printf ( "Rotating phase plane is used.\n"); 

printf ("Real tine FFT function is added .. .\n" ) ; 

printf ( "Energy monitoring available via A H in text node"); 

printf (" A L gives element pinning, NOT A P. . . \n\n\n" ) ; 

printf ("\nDo you want to use a file for initial conditions (Y/N)? "); 

scanf ( "%s" , answer ) ; 
sample__counter=0 ; 

if ( ( answer [ 0 ]“89 ) | | ( answer ( 0 ]==121 ) ) ( /♦ FILE INPUT FOR STARTUP */ 

printf ("Do you want to use old format data file? " ) ; 
scanf ( "%s" , answer ) ; 

printf ( "\nEnter name of file to be read: "); 

scanf ( "ts" , f i lename ) ; 
if ( ( f q=f open( f ilename , "r" ) ) !=NULL) ( 
mean_gamma=0 ; 

f scanf ( f q, "*d\n" , &no_pendulums) ; 
f scanf ( fq,"*lf\n" , talpha ) ; 
f scanf ( fq,"%lf\n" ,6beta) ; 
if ( ( answer ( 0 ]==89 ) | | ( answer ( 0 ]“121 ) ) ( 
f scanf ( f q, "%lf \n" , &mean_gamma ) ; 

DOFOR( i ,no_pendulums) ( 
gamma [ i ]=mean.gamma; 
oraegaO [ i ]=1 ; 

) 

) 

f scanf ( f q , 1 f \n" , 6eta ) ; 
f scanf ( fq, "%lf \n" , Iomega) ; 

DOFOR( i ,no_pendulums) ( 

"• if (( answer ( 0 ] 1=89 )&&( answer [ 0 ] !»121 ) ) { 

f scanf ( fq, "%lf %lf %lf Uf\n'\ 

fcomegaO [ i ) , &gamma [ i ] , ^coordinate ( i ] , &momentum[ i ] ) ; 
mean_gamma=raean_gamma+gamraa[ i ] ; 

) 

else f scanf ( fq, "%lf \lf \n" , &coordinate[ i ] , &momentum[ i ] ) ; 
pinned_elements[ i ]=0; 

) 

if ( ( answer! 0 ] !=89 )&&( answer [ 0 ] !=121 ) ) 
mean_gamma=mean_gamma/no_pendulums; 
f close ( fq ) ; 

) 

else printf ( "Can' t open file requested."); 

) /* if ( ( ans . . . V 
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else < /* USER STARTUP V 

printf ( M \nEnter number of pendulums to use: " ) ? 

scanf ( M td M , &no_pendulums ) ; 

printf ( M \nEnter mode amplitude: " ); 

6canf ( "%lf " , &mode_amp) ; 

printf ( "\nEnter modulation amplitude: ")? 

scanf ("*lf " , *max_amp) ; 

printf ( "\nEnter nonlinear coefficient alpha (+/- 1 ONLY): "); 

scanf ( "%lf 91 , fcalpha ) ? 

DOFOR (k,no_pendulums) ( 

if (alphaco) coordinate! k ]~pow( (-1 ) , k) * (mode_amp 
+max_amp*sin( 2*PI*k/no_pendulums) ) ; 
else coordinated ]*mode_amp+max_amp*sin( 2*PI*k/no_pendulums) 
momentum[k ]=0 ; 
pinned_eleraents[k]»0 ; 

) /* DOFOR */ 

printf ( "\nEnter coupling cpefficient gamma: ” ); 

scanf ( "%lf M , &mean_gamma) ; 

DOFOR ( i , no_pendulums ) { 
gamma [ i ] *=mean_gamma ? 
omegaO[ i ]*1 ; 

) 

printf ( "\nEnter drive amplitude eta: " ) ? 

6canf ( "%lf " , fceta ) ; 

printf ( "\nEnter drive frequency omega: "); 

scanf ( "%lf " , Iomega ) ; 

printf ( "\nEnter dissipation constant beta: "); 

6canf ("%lf",&beta) ; 

) /* else */ 

printf ( "\nEnter time constant (multiple of period/200): ” ); 
scanf ( "%lf " ,&time_int ) ; 
period=2*PI/omega ; 
time_int=time_int*period/200; 

DOFOR ( i , 15 ) flags! i ) =0 ; 

DOFOR ( i , 4096 ) spectrum! i ]=0 ; 
text_f lag=l ; 
mean_omegaO=l . 0 ; 
spectrum_counter=0 ; 

DINC-.02; /* CHARGE THIS TO CHANGE SCALE OF DISPLAY */ 

PINC=2; /* CHANGE THIS TO CHANGE SCALE OF PHASE PLOT */ 

) /* USER_INIT */ 

/**** GRAPHICS DISPLAY ROUTINE ****/ 

display_graphics( ) ( 
int c,i, j,k,l,m,n; 
if (no_pendulums<=40 ) { 
f irst_element=0 ; 

DOFOR ( k , no_pendulums ) ( 
n=0; 

if ( ( stabil i ty_element==k ) && ( stabi li ty_f lag==l ) ) n=l ; 
if ( (pinned_elements[k]“l )&&(disp_color==l ) ) n=2; 
if ( ( pinned_elements [k ] e =l )&( (disp_color=-2 ) ) n=-l ; 
l*old_coord inate [ k ] /DISPLAY_INCREMENT ; 

_setcolor(0) ; 

_setpixel ( ( 160-5*MIDDLE_ELEMENT + 5*k) ,( 100+1 )) ; 
l=coordinatetk]/DISPLAY_INCREMENT; 

_setcolor (disp_color+n ) ; 

_setpixel( ( 160-5* MI DDLE..ELEMENT + 5*k) , ( 100+1 )) ; 
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) /* DOFOR */ 

) /* IF V 

else ( 

1«( no_pendulums~40 )/2 ; 
f irst_element=l ; 

DOFOR( k , 40 ) { 

m=old_coordinate [ 1+k ] /DISPLAY_I NCREMEHT ; 
n=0; 

if ( ( stability _element==( 1+k) ) && ( stabi lity_f lag==l ) ) n=l ; 
if ( (pinned_elements[k]«=l ) ( disp_color==l ) ) n=2; 
if ( (pinned_elements[k]“l )fcfc(disp_color==2) ) n*-l ; 
_setcolor(0) ; 

_setpixel( (60 + 5*k) , ( 100+m) ) ; 
incoordinate [ 1+k ] /DISPLAY_INCREMENT; 

_setcolor ( disp_color+n ) ; 

_setpixel( (60 + 5*k) , ( 100+in) ) ; 

) /* DOFOR */ 

) /* ELSE */ 

) /* DISPLAY_GRAPHICS */ 

/**** TEXT DISPLAY ROUTINE ***♦/ 



display_text( ) ( 

char message [ 80 ] ; 
int c, j ,k, 1 ,m; 

_setvideomode(_DEFAULTMODE) ; 
printf ( "Time is : %lf " ,model_time) ; 
printf (" System parameters are: \n" ); 
printf (" Gamma XI f " ,mean_gamma ) ; 

printf (" Eta *lf",eta); 

printf ( " Omega X 1 f\n" , omega ) ; 

printf (" Beta %lf",beta); 

printf (" Alpha %lf\n",alpha); 

printf ( "\nThere are %d elements in the system" , no_pendulums ) ; 

printf ( "\n\nPress any key to continue..."); 

c*getchar ( ) ; 

while ( kbhit ( )-=0 ) ? 

printf ( "Element Position Velocity Element Position 

DOFOR ( j , 20 ) { 

printf ( " %d %lf %lf %d %lf 

j , coordinate [ f i rst_element + j ) , momentum! f irst_element+ j ) , 

( j+20) , coordinate! f irst_element+ j+20 ) , 
momentum! f irst_element+ j+20 ) ) ; 

) /* DOFOR */ 

printf ( "\n\nPress A G for graphics, A Q to quit." ); 
c=getchar ( ) ; 
while(kbhit( )«0) ; 
while (kbhit( )==0) ; 

) /* display_text */ 



Veloci 

% 



/**★* PAUSE ROUTINE. THIS GIVES USER CHOICE OF SAVING CURRENT STATE 

TO A DATA FILE. ADDITIONALLY, THE TIME STEP CAN BE VARIED DURING 
A PAUSE. *+**/ 



pjrocess_pause( ) { 

int c, i , j ,k, 1 ,mode ; 

FILE *fr; 

char filename! 30 ) , messagefBOJ; 
_clearscreen ( _GCLEARSCREEN ) ? 
_setvideomode (_DEFAULTMODE ) ; 
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printf ("Do you wish to save this state? *' ) ; 
if ( (c=getche( ) )==121 ) ( 

printf ( "\nEnter name of file to be written: " ) ; 
scanf ("%s", filename) ; 
if ( (fr-fopen( filename,*^" ) ) 1-HULL) { 
f printf ( f r , '*%d\n" , no_pendulums ) ; 
f printf ( f r ,"%lf\n" , alpha ) ; 
f printf (fr , "tlfXn" ,beta) ; 
f printf ( fr , "%lf\n" , eta ) ; 
f printf ( f r , "%lf\n" , omega) ; 

DOFOR(i ,stability_ element) 

f printf ( fr , "%lf %lf %lf %lf\n", 

omegaO[ i] ,gamma[ i ) ,old_coordinate[ i ) , momentum! i ] ) ? 

DOFOR( i , ( no_pendulums-stabili ty_element ) ) 

f printf ( fr , "%lf %lf %lf %lf \n" ,omegaO( i+stability_element ] , 
gamma [ i+stability_element ] ^coordinate! i+stability_ element ] , 
momentum! i+stability_element ] ) • 
fclose(f r) ; 

} /* if ( ) */ 

else printf ( "Failed to open %s\n" , filename) ; 

) /* if (()) */ 
time_int=time_int*200/period ; 

printf ( "\nEnter new time multiple (old multiple is %lf): " , time_int) ? 

scanf ( "Ilf” , &time_int) ? 
time_int«time_int*per iod/200 ; 

if (phase_f lag==0) _setvideomode(_MRES4C0L0R) ; 
else _setvideomode (_VRES16C0L0R ) ; 

) /* process_pause */ 

/♦*** START DUMP OF INDIVIDUAL ELEMENT TIME SERIES TO DATA FILE *♦**/ 

start_f ile_dump( ) { 

_clearscreen(_GCLEARSCREEH ) ; 

_setvideomode (_DEFAULTMODE) ; 

printf ( "Enter number of element to be monitored: " ) ; 

scanf ( "%d" , &chosen_element ) ; 

if ( chosen_ element> ( no_pendulums-l ) ) 

printf ("Out of range. No file dump" ) ; 
else f ile_dump_f lag=l ; 

_clearscreen(_ GCLEARSCREEN ) ? 
if ( phase_f lag=*0 ) 

_setvideomode(_MRES4C0L0R) ; 
else 

_ setvideomode (_ VRES16C0L0R) ; 
counterl=0 ; 

) /* start_f ile_dump */ 

/**** COMPLETE DUMP OF INDIVIDUAL ELEMENT TIME SERIES TO DATA FILE ***♦/ 

stop_f ile_dump( ) { 
char filename! 30) ; 
int c; 

FILE *fs; 

_clearscreen(_GCLEARSCREEN) ; 

_setvideomode ( _DEFAULTMODE ) ? 

printf ( "\nEnter name of file to be written: "); 

scanf ("%s" , filename) ; 

if ( ( fs-fopen( filename , "w" ) ) ! -NULL) { 

• .* DOFOR(c, 8000) fprintf ( f s , "%lf\n" , time_series_disp[ c ] ) ? 

_clearscreen(_GCLEARSCREEN) : 
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counter 1=0 ; 
i f ( phase_f lag==0 ) 

_setvideomode(_ MRES4C0L0R) ; 
else _setvideomode(_VRES16C0L0R) ; 

) /* if ( ( ) } V 
) /* stop_f ile_dump */ 

/**** CHECK FOR SYSTEM STABILITY. THIS IS A USER AID ONLY, AND DOES NOT 
PERFORM A RIGOROUS STABILITY CHECK. IT SIMPLY CHECKS TO SEE IF 
ALL OF THE LAST TWENTY PEAK AMPLITUDES OF THE CHOSEN ELEMENTS ARE 
WITHIN ONE PERCENT OF EACH OTHER ♦***/ 

monitor_stability ( ) { 

if (stability_f lag==0) { 
stability_f lag=l ? 
disp_color=2 ; 

_setvideomode(_DEFAULTMODE) ; 

printf ( "Enter number of element to be monitored: n ) ; 
scanf ( "%d" , &stability_element ) ; 
stability_restart (); 

) 

else ( 

stability_f lag=0; 
disp_color=l ; 

) 

if (phase_f lag==0) 

__ setvideomode (_MRES4COLOR) ; 
else 

_setvideomode(_VRES16COLOR) ; 



stability_check( ) { 
int i , k; 
k=0 ; 

DOFOR( i , 19 ) ( 

stable_f lag=l ; 

if ( ( peak__record [ i + 1 ] > ( peak_record ( i )*{ t ♦STAnil.ITY_JNCRF.MrNT) ) ) | 

( penk_record [ i-» 1 ] < ( peak_record [ i ) * ( 1 - STAB1 L1TY__1NCREMENT )))) { 

stable_f lag=0 ; 

break; 

) /* if V 

) /* DOFOR */ 

if (stable_f lag==l ) ( 

disp_color=l; 

) 



/**** pin ANY INDIVIDUAL ELEMENT AT ZERO AMPLITUDE. *♦**/ 

pin_elements( ) { 
int ans,i, check; 

_se tvideomode ( _DEFAULTMODE ) ; 
check=0 ; 

printf ( "\nEnter number of element to be pinned: " ); 

scanf ( "Id" , &ans ) ; 

if ( ( ans<0 ) j ( ans> ( no_pendulums-l ) ) ) ( 
check=l ? 

) 

if(check==0) ( 

if (pinned_elements[ans )==0) < 
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pinned_elements[ ans ]=1 ; 
coordinate[ ans )=0; 
momentum [ ans ] «0 ; 

) 

else pinned_elements[ans]*=0; 
) /* if V 
if (phase_f lag*«=0) 
_6etvideomode(_MRES4COLOR) ; 
else 

_setvideomode(_VRES16COLOR) ; 



/**♦* KICK any individual element any desired amount in amplitude **♦*/ 

kick_elements( ) ( 
int ans , check; 
double amount; 

_setvideomode(_DEFAULTMODE) ; 
check=0; 

printf ( *'\nEnter number of element to kick: "); 

scanf ( "td” , tans) ; 

if ( (ans<0) | (ans>(no_pendulums-l) ) ) ( 
check=l ; 

) 

if(check=*0) ( 

printf ( M \n Enter amount to kick: M ); 

scanf ( "tlf •• , fc amount) ; 

coordinate (ans Incoordinate ( ans] +amount ; 

) /* if V 

if ( phase_f lag”0 ) 

_setvideomode(_MRES4 COLOR) ; 
else 

_setvideomode(_VRES16COLOR) ; 



stability_restart( ) { 
int i; 

disp_color=2 ; 
stablest lag=0 ? 

DOFOR( i , 20 ) peak_record[i)-0; 
peak_record[ 15 ]=10 ; 



/**** INITIALIZE THE PHASE PLANE PLOTTING SYSTEM. THIS ROUTINE ONTAINS THE 

USER CHOICES FOR ELEMENTS TO MONITOR AND DRAWS THE AXES OF THE PLOT. ***♦/ 

init_phasplot( ) { 
int c, i , j ,k, 1 ? 

char ans( 5 ] ,message[ 40 ] , txt[ 3 ] ; 
float dummy? 

spectrum_f lag=0; 

_setvideomode(_DEFAULTMODE) ; 

D0F0R(i,5) phase_elements[ i ]=0; 

printf ( M \nDo you want Poincare sections? M ) ; 

if ( (c=getch( ) )-=121 ) phase_f lag=2 ? 

else phase_f lag=l ; 

printf ( M \nWhich elements do you wish to monitor (999 to finish): ")? 

i— 1; 

k=0? 
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while ( (i++! =999) && (k<5) ) { 

scanf ( "%d" , &i ) ; 
phase_elements [ k++ ]=i + 1 ; 

) /* while */ 

DOFOR(i,5) ( 

if ( phase_elements[ i ]==1000 ) { 
phase_elements[ i ]=0; 

) 

) 

_setvideomode(_VRES 16 COLOR) ; 

_clearscreen (_GCLEARSCREEN ) ; 
graphics_f lag=0; 

_settextcolor(7) ; 
colors [0]*=2; 
colors [ 1 ]=3 ; 
colors [ 2 ]=9 ? 
colorsf 3 ) = 14 ; 
colors[ 4 j=13; 

DOFOR( i , 120 ) ( 

_setpixel( 318, 4*i) ; 

) 

DOFOR( i , 160 ) { 

_setpixel ( 4*i , 238 ) ; 

) 

D0F0R(i,6) { 

_setpixel ( 319 , lO*SCREEN_CORRECTION_FACTOR+l 
+ ( 80/SCREEN_CORRECTION_FACTOR) * ( i+1 ) ) ; 

) 

D0F0R ( i , 8 ) ( 

_setpixel(80*(i+l) ,239) ? 

) 

_inoveto( 480,450)? 

DOFOR( i , 5 ) { 

if ( (phase_elements[ i ) ) !=0) ( 
c=phase_elements[ i )-l ? 

/* itoa(c, txt, 10) ; 

_outtext(txt) ; */ 
printf("%d " ,c); 

_setcolor( colors [ i ]) ? 

DOFOR( j , 4 ) { 

DOFOR(k , 8 ) ( 

_setpixel ( 550+10* i+k , 10+ j ) ; 

) 

) 

) 

) /* DO FOR */ 

printf ( n LAT2X4 86 . C M ) ? 

print f ( M \nGamrna : %lf Eta: %lf Beta: %lf " , mean_gamma , eta, beta) 

printf ( "\nOmega: %lf Alpha: %lf Tine interval: Ilf** , 

onega , alpha ,time_int ) ; 
dummy=l/( time_int*onega ) ? 
phase_int=dummy/l ; 
phase_counter=0 ; 



/**** PHASE PLANE PLOTTING ROUTINE **♦*/ 

phasplot() { 
int i,j,k,l; 

double speed, position, f ast_coordinate, f ast_monentum; 



203 



if ( phasc_f lag==l ) ( 

DOFOR (1,5) | 

if ( ( j-phnr;o_e lemon ts[ i J-l ) !--l ) { 

speed-PIIASE_INCREMENT*momentum( j ) ; 

_setcol or ( colors [ i ] ) ; 

posit ion=PHASE_INCREMENT*coord i nato [ j ) ? 
f ast_coordinate=( position*cos( omega * model _ti me ) 
-speed*sin(omega*model_time) /omega ) 

/ ( 3*SCREEN_C0RRECTI0N_FACT0R) ; 
f as t_momentum=( posit ion*s in ( omega *model_time) 

+speed*cos (omega*roodel_time) /omega )/4 ; 

_setpixel( (320+320*f ast_coordinate ) , ( 240+240*f ast_momentum) ) ; 

) /* if ( ( ) ) */ 

) /* DOFOR */ 

) /* if V 

else if (phase_f lag==2) ( 
if (phase_counter“phase_int) { 
phase_counter*0 ; 

DOFOR ( i , 5 ) { 

if ( ( j=phase_elements [ i ]-l) !*-l) ( 

speed=PHASE_INCREMENT*momentum [ j ) ; 

_setcolor ( colors [ i ] ) ; 

position=PHASE_INCREMENT*coordinate[ j ] ; 
f ast_coordinate=(position*cos(omega*model_time) 
-speed*sin(omega*model_time) /omega ) ; 
f ast_momenturo=( position*sin( omega*model_time ) 

+speed*cos ( omega*model_ time ) /omega ) /4 ; 

_ setpixel ( ( 320+320* f as t_coordinate) , ( 240+24 0*f ast_jnomentum) ) ; 

) /* if ( ( ) ) V 
) /* DOFOR */ 

) /* if V 

) /* else if ( ) */ 

) /* phasplot */ 

/* FFT ROUTINE */ 

compute_spectrum( ) { /* Uses algorithm from p. 411 of Numerical Recipes in C */ 
int n, inmax ,m, j , istep, i ,nn, temp; 
double wtemp,wr ,wpr ,wpi ,wi , theta , tempr , tempi ; 
nn=2048 ; 
n*nn«l ; 

j-i; 

DOFOR ( i , 204 8 ) { 

nn»(i&1024)/1024; 

nn=nn+(i&512 )/256; 

nn=nn+( i&256 ) /64 ; 

nn=nn+( i&128 )/16; 

nn=nn+( i&64 )/4 ; 

nn=nn+(i&32) ; 

nn=nn+( i£16) *4 ; 

nn=nn+( i&8) *16; 

nn=nn+( i&4 )*64 ; 

nn=nn+(i&2 )*256; 

nn=nn+( i&l)*1024; 

temp2[ 2*nn J=spectrum[ 2 * i J ; 

temp2[ 2*nn+l ]=spectrum[ 2*i+l ] ; 

) 

DOFOR( i , 4 096 ) spectrum[i ]=temp2[i ] ; 
mmax=2 ; 

while (n>mmax) ( 



204 



istep^2*mmax ; 
theta=2*PI/mmax ; 
wtemp=sin ( 0 . 5* theta ) ; 
wpr=-2 . 0*wtemp*wtemp; 
wpi^sin(theta) ; 
wr=l . 0 ; 
wi«0 . 0 ; 

f or(m*0 ;m<(mmax-l ) ;m+»=2 ) ( 
for ( i=m? i<=n; i+*istep) { 
j=i+mmax; 

tempr=wr ♦spectrum [ j ]-wi*spectrum[ j + 1 ] ; 
tempi=wr*spectrum[ j+1 ]+wi*spectrum[ j j ; 
spectrum[ j ]=spectrura[ i ]-tempr; 
spectrum[ j+1 ]=spectrum[ i+1 ]-tempi ; 
spectrum [ i ]+=tempr ; 
spectrum[ i + 1 ]+=tempi ; 

) 

wr=(wtemp=wr ) *wpr-wi *wpi+wr ; 
wi=wi*wpr+wtemp*wpi+wi ; 

) 

nunax-istep? 

) 

) 

/**** SPECTRUM PLOTTING ROUTINE *♦**/ 

plot_spectrum( ) { 
int i , j ; 

double freq , magnitude , mag , max; 
int c; 

FILE *f p ; 

if (dump_spectrum_f lag==l ) ( 

printf ( “\nDumping spectrum now" ) ; 
f p-f open ( “spectrum. out” , M w” ) ; 
for(i=l ;i<204 8 ?i++) ( 

if ( time_int !«0) f req=i/( 2048*time_int ) ; 
else printf ( ”Time_int was zero!")? 
magnitude=sqrt (spectrum[ ( 2*i) ]*spectrum[ ( 2*i ) ] 

+spectrum[ ( 2*i + l ) ] *spectrum[ (2* i+1)]); 
f printf (fp, "%lf %lf \n” , freq, magnitude) ; 

) 

f close( f p) ; 
durop_spec trum_f 1 ag=0 ; 

) 

_clearscreen(_GCLEARSCREEN) ; 
max=0; 

for (i=l ; i<2049;i++) { /* this loop calculates, the max spectral component*/ 
magnitude=sqrt (spectrumf ( 2*( i ) ) ]*spectrum[ ( 2* (i ) ) ] 

+spectrum[ ( 2* ( i )+l ) ] *spectrum[ (2*(i)+l)j); 
if (magnitude>max) max=magnltude; 

) 

_setcolor ( 3 ) ; 

DOFOR(i,640) _set pixel ( i , 460 ) ; 

DOFOR ( i , 1 1 ) _setpixel(51*i,461 ) ,_setpixel ( 51 *i ,462) ; 
printf ( “Spectrum for element td LATTIC2X . C“ , spectruro_element ) ; 

printf ( “\nAlpha : %lf Beta: %lf Gamma: %lf Eta: tlf Omega: %lf“, 
alpha, be ta ,mean_gamma , eta, omega) ; 
printf ( “\nMax amp = %lf M ,max); 

freq=512/( 204 8* time_int ) ; 

printf (“Max freq shown: %lf Time interval: tlf “ , f req , tirae_int ) ; 
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_setcolor ( 1 ) ; 

DOFOR(i,512) ( 

j-o; 

©agnitude=sqrt( spectrum [ 2*i )*spectrum[ 2*i ] 
♦spectrum [ 2*i + l ] * spectrum [ 2*i+l J ) ? 
jnag=magnitude*4 00/max; 
while( j++<mag) ( 

_setpixel( i , 460- j) ; 

} /* while */ 

) /* DOFOR */ 



/**** INITIALIZES WATERFALL DISPLAY MODE, DRAWING GRID **♦*/ 

init_waterfall ( ) ( 
int c,i,j,k,l,m; 
char ans[5]; 

_setvideomode (_VRES1 6COLOR ) ; 

_setcolor ( 1 ) ; 

DOFOR ( i ,600 ) _setpixel ( 20+i ,20); 

DOFOR ( i , 450 ) _setpixel ( 20 , i+20 ) ? 

DOFOR( i , 61 ) ( 

if ( ( i==10 ) | (i==20) | (i==32) | (i==42) | (i==52) ) _setcolor( 2 ) ; 
DOFOR ( j , 113 ) { 

_setpixel (20+10*i , 20+4* j ) ; 

) 

_setcolor(l); 

) 

_setcolor ( 2 ) ; 

DOFOR ( i , 1 1 ) ( 

DOFOR( j , 160 ) ( 

_setpixel(20+4* j, 16 + 40* (i + 1) ) ? 

) 

) 

waterf a 1 l_count er=0 ; 
col or. count er=2 ? 



/*♦** WATERFALL DISPLAY ROUTINE ****/ 

disp_waterfall ( ) ( 
int c,i, j,k,l; 

if (waterf all_counter++<110) { 
if (color_counter++==14 ) { 
color_counter=2 ; 

_setcolor(color_counter) ; 

) 

else ( 

_setcolor(color_counter ) ; 

) 

DOFOR ( i ,no_pendulums) ( 
if (waterf all_counter<55 ) ( 

_setpixel ( 20+2*i , 24+8*waterf al l_counter-7*wat_inc* (coordinate! i ) ) ) 

) 

else ( 

_setpixel ( 340+2*i , 24 +8* (waterf all_counter-55 ) -7 *wat_inc*( coordinate [ i ) ) ) 
) 

) 

) 

else ( 
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